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Abstract 

The  Kolmogorov  Arnold  and  Moser  (KAM)  theorem  states  that  a  lightly  per¬ 
turbed  Hamiltonian  system  will  have  solutions  which  lie  on  a  torus.  The  trajectories 
of  Earth  orbiting  satellites  have  been  shown  to  lie  on  KAM  tori  with  three  basis  fre¬ 
quencies.  These  basis  frequencies  are  determined  by  fitting  second  order  polynomials 
to  data  from  Two-Line  Element  sets  (TLEs)  using  a  least  squares  technique.  The 
first  order  coefficients  are  used  as  the  torus  basis  frequencies  while  the  second  order 
terms  are  used  to  account  for  perturbations  to  the  satellite’s  orbit  such  as  air  drag. 
Four  cases  are  attempted  using  the  Hubble  Space  Telescope  and  three  rocket  bodies 
as  test  subjects.  A  KAM  torus  with  the  desired  basis  frequencies  is  constructed  and 
used  to  predict  satellite  position.  This  position  prediction  is  compared  to  the  position 
obtained  from  TLEs  using  Simplified  General  Perturbations  4  (SGP4)  algorithms. 
Analysis  of  the  torus  position  error  shows  that  the  torus  construction  algorithm  does 
not  fully  characterize  the  contribution  of  the  smaller  basis  frequencies  to  the  orbital 
trajectory. 
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KAM  Torus  Frequency  Generation  from  Two-Line 

Element  Sets 

I.  Introduction 

1 . 1  Motivation 

Since  the  launch  of  Sputnik  by  the  Soviet  Union  in  1957,  there  has  been  a 
steady  increase  in  the  importance  of  space-based  assets  to  everyday  life.  Today, 
everything  from  banking  to  car  navigation  to  watching  sports  on  television  is  either 
directly  tied  to,  or  influenced  by  capability  provided  by  satellites.  The  increasing 
importance  of  Space  has  not  been  limited  to  the  civilian  world.  In  the  past  decades, 
the  U.S.  military  has  come  to  increasingly  rely  on  capabilities  provided  by  space 
assets  in  everyday  operations.  These  assets  provide  vital  services  such  as  communi¬ 
cations,  imagery  and  signals  intelligence,  and  global  positioning  data  to  U.S.  Military 
personnel  and  our  allies  around  the  world. 

In  order  to  communicate  with  these  on-orbit  assets,  fixed  ground  sites  and 
mobile  users  must  know  where  a  satellite  will  be  at  a  given  point  in  time.  Because  of 
this  requirement,  it  is  vital  that  we  have  an  accurate  picture  of  where  all  space  objects 
are  located  at  the  current  time,  and  the  capability  to  rapidly  and  accurately  predict 
where  they  will  be  in  the  future.  In  addition  to  our  own  satellites,  the  near-earth 
environment  is  filled  with  thousands  of  operational  commercial  and  foreign  owned 
satellites  as  well  as  many  more  thousands  of  pieces  of  space  junk  including  dead 
satellites,  spent  rocket  bodies,  debris  from  orbital  collisions  and  trash  from  past 
manned  space  missions.  These  objects  present  a  hazard  to  operational  satellites, 
especially  when  their  precise  positions  are  unknown. 

The  Air  Force’s  current  method  of  orbit  propagation  uses  the  Simplified- 
General  Perturbations  4  (SGP4)  model  which  dates  back  to  the  1970s.  This  model 
has  been  adequate  but  can  only  accurately  predict  a  satellite’s  position  for  a  few 
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days  before  an  update  is  necessary.  If  more  accurate  predictions  are  needed,  for 
collision  avoidance  for  example,  numerical  integration  is  used.  While  advances  in 
computational  technology  have  greatly  reduced  the  time  necessary  to  perform  these 
integrations,  they  still  take  hours  to  complete.  In  addition,  to  predict  position  one 
week  in  advance,  the  orbit  must  be  integrated  through  all  the  intermediate  time 
steps  to  produce  the  one  needed  predicted  location.  A  new  method  with  greater 
accuracy  and  less  computational  and  observational  cost  would  be  beneficial. 

1 . 2  Approach 

This  effort  will  study  the  feasibility  of  converting  the  current  method  of  orbit 
propagation  using  SGP4  to  one  involving  a  KAM  torus.  Many  of  the  methods  used 
in  the  construction  of  the  tori  are  based  on  those  developed  and  demonstrated  in  the 
past  by  Wiesel  and  his  past  students.  This  approach  may  provide  increased  orbital 
prediction  accuracy  at  a  significantly  lower  computational  and  observational  cost. 

1.3  Problem  Statement 

This  work  will  answer  the  question  of  whether  distinct  basis  frequencies  can 
be  extracted  from  Two-Line-Element-Set  (TLE)  data,  whether  KAM  tori  with  those 
basis  frequencies  can  be  constructed  from  initial  position  and  velocity  obtained  from 
TLEs  and  SGP4,  and  whether  these  tori  can  accurately  predict  future  positions  of 
an  Earth  orbiting  satellite  at  a  level  of  accuracy  equal  or  greater  than  that  of  SGP4. 

1.4  Results 

KAM  torus  frequencies  were  identified  for  two  of  four  test  cases  showing  that, 
for  certain  satellite  types  in  certain  orbits,  it  is  possible  to  extract  frequencies  from 
purely  observational  data.  In  the  two  unsuccessful  cases,  it  is  believed  that  a  com¬ 
bination  of  variations  in  air  drag  along  with  inaccuracies  in  the  TLE  data  were  the 
cause  of  the  poor  curve-fits  to  TLE  data.  KAM  torus  construction  was  attempted 
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for  the  remaining  two  test  cases.  Of  these,  the  construction  process  was  unsuccessful 
for  one  case,  most  likely  clue  to  the  very  small  eccentricity  of  the  orbit.  In  the  final 
case,  a  KAM  torus  was  constructed  and  its  predicted  position  was  compared  to  the 
position  predicted  at  each  TLE  epoch  using  SGP4  algorithms.  This  position  com¬ 
parison  showed  promising  results  in  that  the  position  error  of  the  torus  prediction 
showed  very  little  linear  growth.  The  error  did,  however,  show  significant  periodic 
oscillation.  These  oscillations  were  shown  to  occur  at  the  two  smallest  torus  basis 
frequencies.  It  is  believed  that  this  is  caused  by  an  error  in  the  torus  construction 
process,  namely  the  calculation  of  Fourier  coefficients. 
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II.  Background 


2.1  Space  Surveillance  Network 


The  Space  Surveillance  Network  (SSN)  is  a  network  of  29  sites  located  around 
the  world.  These  sites  detect  and  track  man-made  objects  orbiting  the  Earth  in¬ 
cluding  operational  and  non-operational  satellites,  spent  rocket  bodies,  debris,  and 
fragments  [28].  The  SSN  operates  a  variety  of  sensors  to  accomplish  its  mission  in¬ 
cluding  phased-array  and  conventional  radars,  electro-optical  sensors,  the  Mid  course 
Space  Experiment  (MSX)-  a  Low  Earth  Orbit  (LEO)  satellite  with  a  payload  of  sen¬ 
sors  spanning  UV  to  very-long-wave  IR  -  and  ground-based  electro-optical  deep  space 
surveillance  sites  which  provide  tracking  of  deep  space  objects,  including  geostation¬ 
ary  satellites,  orbiting  above  22,500  miles.  Combined,  these  sensors  are  responsible 
for  300,000  to  400,000  observations  each  day  [28].  Figure  1  shows  how  these  sensors 
are  spread  over  the  Earth  to  provide  maximum  coverage. 


Space  Surveillance  Network 

Worldwide  Network  of  20  Optical  and  Radar  (Mechanical  &  Phased  Array)  Sensor  Sites 
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Figure  1:  SSN  Sites  [27] 


The  SSN  monitors  space  objects  using  a  predictive  technique  rather  than  a  con¬ 
tinuous  approach.  This  means  because  of  the  number  of  objects  being  tracked  and 
the  limited  capability  of  instruments,  objects  are  not  tracked  in  real  time.  Rather, 
their  position  is  acquired  only  periodically.  This  position  is  then  propagated  for¬ 
ward  in  time  using  the  SGP4  dynamics  model.  At  some  later  time,  SSN  instruments 
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then  look  for  the  object  at  the  position  it  was  predicted  to  be  and  records  its  ac¬ 
tual  position,  which  is  hopefully  close  to  the  prediction.  If  a  satellite  has  made  an 
unexpected  maneuver  during  the  time  between  observations,  or  if  the  object  has 
experienced  greater  than  expected  perturbing  forces  such  as  more  air  drag  due  to  a 
change  in  the  Earth’s  outer  atmosphere,  the  process  of  re-acquiring  the  object  can 
become  time-consuming.  A  more  detailed  explanation  of  orbit  determination,  and 
of  SGP4,  can  be  found  in  Section  2.5. 


2.2  Orbital  Debris 

Since  the  beginning  of  the  space  age,  the  number  of  objects  in  orbit  has  steadily 
increased  to  the  current  number  of  over  15,000  objects  10cm  in  diameter  or  greater 
being  tracked  by  the  SSN.  Figure  2  provides  a  snapshot  of  the  growth  of  near-earth 
objects  being  tracked  by  the  SSN.  The  large  jumps  in  total  number  of  tracked  objects 
(in  2007  and  2009  for  example)  are  results  of  satellites  breaking  up  due  to  explosions 
or  collisions,  both  intentional  and  unintentional. 


Figure  2:  Plot  of  the  Amount  of  Debris  >10cm  in  Diameter  Being  Tracked  by  the 
SSN  [23] 
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Table  1  shows  the  top  10  such  break-ups,  in  terms  of  debris  created.  Ap¬ 
proximately  one  third  of  all  objects  currently  being  tracked  are  a  result  of  these 
break-ups.  These  collisions  and  explosions  result  in  debris  trains  that  can  rapidly 
circle  the  earth,  endangering  any  near-by  satellites  and  effectively  prohibiting  any 
operational  satellite  from  occupying  or  crossing  that  orbital  plane.  This  is  demon¬ 
strated  by  Figure  3  which  shows  the  debris  train  as  of  December  2007  created  by  the 
January  2007  Chinese  Anti-Satellite  (ASAT)  test  which  destroyed  the  FENGYUN 
1C  polar-orbiting  weather  satellite. 


Table  1:  Top  10  Breakups  as  of  May  2010  [24] 


Common  Name 

Year  of 
Breakup 

Altitude  of 
Breakup  (km) 

Cataloged 

Debris 

Debris 

in  Orbit 

Cause  of 
Breakup 

Fengyun-lC 

2007 

850 

2841 

2756 

Intentional 

Collision 

Cosmos  2251 

2009 

790 

1267 

1215 

Accidental 

Collision 

STEP  2 

Rocket  Body 

1996 

625 

713 

63 

Accidental 

Explosion 

Iridium  33 

2009 

790 

521 

498 

Accidental 

Collision 

Cosmos  2421 

2008 

410 

509 

18 

Unknown 

SPOT  1 

Rocket  Body 

1986 

805 

492 

33 

Accidental 

Explosion 

OV  2-1/LCS  2 
Rocket  Body 

1965 

740 

473 

36 

Accidental 

Explosion 

Nimbus  4 

Rocket  Body 

1970 

1075 

374 

248 

Accidental 

Explosion 

TES 

Rocket  Body 

2001 

670 

370 

116 

Accidental 

Explosion 

CBERS  1 

Rocket  Body 

2000 

740 

343 

189 

Accidental 

Explosion 

Total:  7903 

Total:  5172 

The  hundreds  to  thousands  of  pieces  of  debris  caused  by  each  break-up  are 
completely  uncontrolled,  and  therefore  present  a  significant  danger  to  operational 
satellites  orbiting  nearby.  As  the  number  of  orbital  objects  increase,  the  probability 
of  future  collisions  increases.  This  phenomenon  has  been  examined  by  Kessler  who 
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Figure  3:  FENYUN  1C  Debris  (red)  from  Chinese  ASAT  Test  and  ISS  Orbit 
(green)  [6] 

predicts  an  increase  in  collisions  and  therefore  an  increase  in  the  amount  of  debris  to 
a  point  when  the  rate  of  new  debris  being  created  exceeds  the  rate  at  which  objects 
enter  the  atmosphere  and  therefore  are  eliminated  [16].  At  that  point,  entire  regions 
of  space  could  become  unusable  for  operational  spacecraft.  While  many,  in  fact  the 
majority,  of  objects  in  orbit  are  uncontrolled,  those  satellites  that  are  most  crucial 
are  those  that  are  currently  in  operation  and  can  therefore,  for  the  most  part,  be 
maneuvered.  A  more  accurate  method  of  predicting  orbital  positions  into  the  future 
will  enhance  the  ability  of  these  satellites  to  avoid  collisions. 

2.3  Collision  Avoidance 

The  importance  of  precise  knowledge  of  a  satellite’s  position  to  collision  avoid¬ 
ance  operations  can  be  seen  from  the  following  example.  Given  a  1km  uncertainty 
in  position  for  two  satellites  with  lm2  cross  sections,  there  is  an  approximately  1  in 
1,000,000  chance  that  the  satellites  will  actually  collide.  Therefore,  if  the  satellites 
were  maneuvered  each  time  a  collision  was  possible  they  would  maneuver  almost 
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1,000,000  times  to  avoid  possible  collisions  before  they  would  maneuver  to  prevent 
an  actual  collision.  The  vast  majority  of  the  time,  in  other  words,  the  satellites  would 
pass  each  other  with  plenty  of  room  to  spare  even  when  a  collision  was  predicted.  On 
the  other  end  of  the  spectrum,  the  Cosmos-Iridium  collision  of  2009  demonstrates 
the  fact  that  predictions  of  a  safe  close  approach  can  be  unreliable.  In  that  instance, 
Satellite  Orbital  Conjunction  Reports  Assessing  Threatening  Encounters  in  Space 
(SOCRATES)  predicted  a  close  approach  of  584  m  between  Iridium  33  and  Cosmos 
2251  only  two  hours  before  they  collided  [13]. 

The  discussions  of  Sections  2.1  through  2.3  demonstrate  that  a  more  accurate 
and  timely  method  of  orbit  determination  would  be  very  beneficial  -  both  in  the 
computational  and  time  savings  realized  in  maintaining  the  orbital  catalog  and, 
more  importantly,  in  providing  data  that  could  be  used  to  more  accurately  conduct 
collision  avoidance  for  U.S.  Military  satellites.  The  recent  and  current  research  by 
Wiesel  and  his  students  is  attempting  to  determine  if  a  method  based  on  KAM 
theorem  could  possibly  meet  this  need. 

2.4  Historical  U.S.  Orbit  Determination 

The  U.S.  Government  interest  in  tracking  space-objects  stems  from  civilian, 
scientific,  and  military  needs.  The  Air  Force  initially  backed  efforts  to  track  and 
catalog  space  objects  in  order  to  have  situational  awareness  in  the  near-earth  envi¬ 
ronment  and  therefore  have  the  ability  to  distinguish  between  a  harmless  orbiting 
satellite  and  a  hostile  incoming  missile.  The  first  effort  to  formally  track  and  catalog 
space  objects  took  place  at  the  National  Space  Surveillance  Control  Center  (NSSCC) 
at  Hanscom  Field  in  Massachusetts.  Satellite  observations  were  taken  at  more  than 
150  sites  around  the  world  using  instruments  such  as  radar,  cameras,  telescopes, 
and  radio  receivers.  The  observational  data  was  all  sent  back  to  the  NSSCC  where 
they  were  first  processed  manually  to  calculate  corrections  to  orbital  elements  before 
being  fed  into  a  computer  to  produce  updated  orbital  data  and  ephemeris.  The  com- 


puter  output  was  sent  back  to  the  observation  sites  in  the  form  of  bulletins  which 
included  data  for  the  next  3-7  days  and  were  used  by  the  sites  to  plan  future  obser¬ 
vation  opportunities.  These  bulletins  were  the  predecessors  to  the  TLEs  currently 
in  use  [12]. 

The  U.S.  Navy  started  development  on  the  Naval  Space  Surveillance  System 
(NAVSPASUR)  in  1958.  This  system  consists  of  a  continuous-wave  multistatic  radar 
interferometer  including  3  transmitters  and  6  receivers  spread  out  along  a  great-circle 
arc  across  the  country  from  San  Diego,  CA  to  Savannah,  GA  and  is  today  commonly 
known  as  “The  Fence” .  This  system  is  unique  in  that  it  detects  all  objects  passing 
through  its  sensing  area  without  any  prior  knowledge  of  the  object’s  orbit.  When 
The  Fence  became  operational  in  1961,  a  computer  required  15  minutes  to  update  a 
single  orbit  [12]. 

The  early  efforts  at  orbital  tracking  by  both  the  Air  Force  and  Navy  used 
a  highly  simplified  dynamics  model  for  orbit  determination.  In  1959  Brouwer  and 
Kozai  published  solutions  for  motion  of  a  satellite  under  the  influence  of  the  Earth’s 
zonal  harmonics  [3]  [17].  These  solutions  did  not  include  air  drag  and  were  con¬ 
sequently  modified  using  various  techniques  concluding  with  Lane  and  Cranford’s 
modification  of  the  Brouwer  solution  in  1969  [18].  Brouwer’s  solution  was  formu¬ 
lated  in  terms  of  Delaunay  variables  and  therefore  contained  the  familiar  problem  of 
small  divisors  with  eccentricity  and  inclination.  This  problem  was  solved  by  Lyddane 
in  1963  who  showed  the  the  Brouwer  solution  could  also  be  expressed  in  terms  of 
Poincare  variables  which  are  the  canonical  variable  counterparts  to  the  equinoctial 
orbital  elements  and  therefore  contain  singularities  only  for  retrograde  equatorial 
orbits  [21].  Brouwer’s  solution  with  Lyddane’s  modifications  was  adopted  by  NAVS¬ 
PASUR  in  1964  and  became  known  as  the  Position  and  Partials  as  functions  of  Time 
(PPT)  model.  The  PPT  model  included  an  approximation  of  air  drag  modeled  as 
influencing  the  mean  motion  as  a  quadratic  function  of  time.  In  1997,  effects  from 
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the  sun  and  moon  were  added  to  the  PPT  model  for  high-altitude  satellites  and  the 
model  became  known  as  PPT3  which  is  currently  in  use  [12]. 

The  Air  Force  also  adopted  a  model  based  on  the  work  of  Brouwer  and  Kozai 
starting  in  1960.  This  model  used  an  air  drag  model  similar  to  PPT,  but  instead 
of  using  Lyddane’s  method  to  avoid  the  small  divisors  problem,  periodic  terms  con¬ 
taining  eccentricity  were  neglected.  This  model  became  known  as  the  Simplified 
General  Perturbations  (SGP)  model  and  became  operational  in  1964.  By  1969,  the 
number  of  satellites  had  increased  to  the  point  of  straining  the  ability  of  current-day 
computers  to  handle  the  numerous  terms  in  the  SGP  model.  This  led  to  the  empir¬ 
ical  atmospheric  density  model  being  replaced  with  an  analytic  one  and  SGP  being 
re-worked  based  on  Lane  and  Cranford’s  work  culminating  in  a  new  model,  SGP4, 
becoming  fully  operational  in  1979  [12]. 


2.5  Current  Orbit  Prediction  and  TLE  Generation 

The  SGP4  model  is  used  in  conjunction  with  TLEs.  A  TLE  is,  as  the  name 
implies,  two  lines  of  data  providing  information  such  as  object  identification,  time  of 
observation,  and  orbital  elements.  Figure  4  shows  an  example  of  a  single  TLE  from 
the  Hubble  Space  Telescope  (HST)  along  with  a  description  of  each  element. 
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Figure  4:  Example  TLE  for  the  HST 


TLEs  for  space  objects  in  Low-Earth  Orbit  (LEO)  are  generated  following 
initial  orbit  determination  from  observations  taken  by  SSN  instruments.  SGP4  is 
used  for  objects  with  periods  less  than  225  minutes,  while  objects  with  orbital  periods 
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greater  than  225  minutes  use  the  Simplified  Deep  Space  Perturbations  4  (SDP4) 
dynamics  model  which  includes  additional  effects  such  as  perturbations  due  to  the 
Sun  and  Moon.  The  SGP4  model  will  be  discussed  here  since  all  orbits  analyzed 
for  this  effort  are  “near-Earth”  orbits.  As  discussed  in  Section  2.4,  the  SGP4  model 
was  originally  based  on  a  theory  of  satellite  motion  perturbed  by  the  Earth’s  zonal 
harmonics  developed  by  Brouwer  [3].  The  model  is  initialized  with  the  following 
parameters,  all  of  which  are  contained  within  a  TLE: 

t0  -  epoch  time 
no  -  mean  motion  at  epoch 
eo  -  eccentricity  at  epoch 
in  -  inclination  at  epoch 

(1) 

lo0  -  argument  of  perigee  at  epoch 
Do  -  right  ascension  of  the  ascending  node  at  epoch 
Mo  -  mean  anomaly  at  epoch 
B*  -  atmospheric  drag  coefficient 

All  of  the  orbital  elements  in  Equation  1,  except  the  mean  motion,  are  mean  orbital 
elements  defined  by  Brouwer  [3].  Brouwer’s  mean  elements  are  based  on  the  first  five 
terms  in  the  Geopotential  (,/]  through  J5)  and  contain  short-period  and  long-period 
oscillations.  Short  period  terms  contain  the  mean  anomaly  in  their  arguments  while 
long-period  terms  contain  multiples  of  the  mean  argument  of  perigee  in  their  argu¬ 
ment.  The  mean  motion  follows  the  convention  developed  by  Kozai  and  includes 
only  short-period  oscillations  [17].  Neglecting  long  period  perturbations  for  mean 
motion  can  be  a  reasonable  assumption  because,  for  many  earth  orbits,  long  period 
effects  are  masked  due  to  the  influence  of  air-drag,  which  decreases  the  orbital  period 
and  therefore  increases  the  mean  motion  over  long  time-scales.  Given  these  initial 
parameters,  the  SGP4  model  can  be  used  to  propagate  the  orbital  elements  forward 
in  time  by  taking  into  account,  again,  the  Earth’s  zonal  harmonics  and  also  atmo- 
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spheric  drag  which  is  modeled  based  on  a  power-law  density  function.  For  a  detailed 
description  of  the  SGP4  equations,  see  Appendix  B  in  [12], 

The  process  of  updating  a  TLE  after  initial  orbit  determination  starts  from 
a  new  set  of  SSN  observations,  each  with  an  associated  time.  A  typical  SSN  radar 
observation  contains  range,  range-rate,  azimuth  and  elevation  while  an  optical  ob¬ 
servation  contains  only  angular  data.  These  observations  then  are  compared  to  the 
predicted  observations  of  the  satellite.  The  predicted  values  are  calculated  by  prop¬ 
agating  the  orbit  forward  in  time  from  the  previous  TLE  using  SGP4.  Then,  the 
predicted  position  and  velocity  of  the  satellite  can  be  determined  at  each  new  obser¬ 
vation  time  and  this  information  can  be  converted  to  predicted  observation  values 
for  the  particular  sensor  which  made  the  actual  observation.  Having  now  obtained 
both  predicted  and  actual  observation  data,  the  goal  is  to  minimize  the  difference 
between  these  two  sets  of  values.  This  is  done  via  a  process  called  differential  cor¬ 
rections.  Let  G  be  an  observation  function  which  expresses  the  observed  values,  z, 
in  terms  of  the  TLE  values,  x 

z  —  G(x)  (2) 

What  is  needed  now  is  a  characterization  of  the  effect  of  changing  a  TLE  value, 
Xi ,  on  the  observation.  This  could  be  accomplished  by  taking  a  derivative  of  the 
observation  function,  G ,  however  is  complex  and  not  easily  differentiable  as  it 

takes  initial  TLE  values,  propagates  them  forward  using  SGP4,  and  then  transforms 
them  to  appropriate  observation  data.  Therefore,  an  approximation  of  the  derivative 
can  be  made  based  on  the  effect  of  a  small  change  in  x,  Ax: 


dG(x )  G(x  +  Ax)  —  G(x) 

dx  Ax 


(3) 
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With  this  approximate  derivative,  new  values  of  x  can  be  calculated  that  minimize 
the  difference  between  the  predicted  and  actual  observations. 

(4) 

X=Xp 

Once  a  new  TLE  is  calculated,  xn,  the  process  is  repeated  replacing  the  predicted 
TLE,  xp  with  the  new  value  until  the  difference  between  the  two  reaches  some  small 
tolerance.  For  a  more  rigorous  discussion  of  differential  corrections  and  initial  orbit 
determination,  see  Wiesel  or  Vallado  [36]  [30].  New  TLEs  are  generally  issued  if 
the  difference  in  predicted  position  between  the  old  and  new  element  sets  is  greater 
than  5km  [14]. 

2.6  Accuracy  of  TLEs  and  SGP4 

The  accuracy  of  a  given  TLE  varies  depending  on  the  number  of  observations 
used  to  generate  the  TLE,  the  accuracy  of  each  observation,  the  current  space  en¬ 
vironment,  the  type  of  orbit,  and  other  factors.  In  general  however,  a  given  TLE  is 
typically  only  valid  for  a  few  days  at  the  most  before  it  needs  to  be  replaced.  Kelso 
compared  the  position  of  GPS  satellites  derived  from  TLEs  and  SGP4  to  precision 
ephemeris  data  and  found  the  TLE/SGP4  position  to  be  accurate  to  only  approxi¬ 
mately  10km  over  a  period  of  15  days  [15].  The  velocity  predictions  from  TLEs  and 
SGP4  are  less  accurate  than  the  position  because  velocity  calculations  rely  on  the 
rates  of  change  of  the  orbital  parameters  and  the  SGP4  model  uses  certain  assump¬ 
tions,  such  as  a  truncated  geopotential  model,  to  determine  these  rates  of  change. 

2.7  Earth  Gravitational  Model  1996 

The  gravitational  model  used  in  this  effort  to  obtain  orbital  position  data 
through  numerical  integration  is  the  Earth  Gravitational  Model  1996  (EGM96)  de¬ 
veloped  by  the  National  Imagery  and  Mapping  Agency  (NIMA),  NASA,  and  the  Ohio 


x, 


Xp  T  za  G(xp) 


dG(x) 
dx 
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State  University.  This  model  was  developed  using  surface  gravity  data  collected  by 
NIMA  as  well  as  satellite  altimetry  and  tracking  data  from  over  20  satellites  and 
consists  of  spherical  harmonic  coefficients  complete  to  order  and  degree  360  [10]. 
Figure  5  shows  a  visual  representation  of  EGM96  in  the  form  of  deviation  from  the 
Earth’s  Geoid  (a  sphere  with  radius  equal  to  the  mean  ocean  surface  of  the  earth). 


Figure  5:  Visual  Representation  of  EGM96  [10] 


The  EGM96  coefficients  are  used  to  construct  the  expression  for  the  Earth’s 
gravitational  potential,  V,  in  spherical  coordinates 

oo  oo  ✓  \  —  n 

V  ( r ,  X,S)  =  -~y]  ^  (  j—  )  P™(cos(<5))  [Cnm  cos  (mA)  +  Snm  sin  (mA)]  (5) 

n=0  m= 0  '  ®  ' 

which  satisfies  Laplace’s  equation: 

VV(r,A,i)  =  0  (6) 

where  r  is  the  scalar  radius,  5  is  the  latitude,  and  A  is  the  east  longitude.  In 
Equation  5,  Cnm  and  Snm  are  the  spherical  harmonic  coefficients  from  EGM96  and 
P™( cos(5))  are  the  associated  Legendre  polynomials  in  cos(<5).  For  a  derivation  of 
Equation  5  see  Wiesel  [35].  For  this  effort,  only  terms  of  order  and  degree  less  than  20 
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in  the  geopotential  were  used  to  allow  computations  to  be  executed  in  a  reasonable 
amount  of  time. 


2.8  Hamiltonian  Dynamics 

Hamiltonian  dynamics  is  a  reformulation  of  classical  dynamics  first  introduced 
in  1833  by  William  Hamilton.  Hamilton’s  methods  can  be  used  to  simplify  a  complex 
dynamical  problem  by  writing  the  equations  of  motion  as  first  order  differential 
equations.  For  a  given  dynamical  system,  the  Lagrangian,  L ,  can  be  written  as  a 
expression  of  the  energy  in  the  system  in  terms  of  the  system  coordinates,  their  time 
derivatives,  and  time 

L(qi,qht)  =  T-V  (7) 

where  T  is  the  system  kinetic  energy  and  V  is  the  system  potential  energy.  L  then 
satisfies  Lagrange’s  equations  of  motion: 

d  f dL\  d L 
dt  \  d<ji )  dq{ 


which  are  second  order  differential  equations.  The  conjugate  momenta  for  the  coor¬ 
dinates  qt  can  the  be  written  as 


Pi  = 


dL 

drji 


(9) 


A  Hamiltonian,  PL,  can  then  be  written  for  the  system 


' H(qi,qi,Pi,t )  =  ^ P&  ~  L  (10) 

i 

The  velocities,  q,  can  be  re-written  in  terms  of  q% ,  pt,  and  t  by  inverting  the  expres¬ 
sions  for  the  momenta  from  Equation  9  and  the  Hamiltonian  can  then  be  written  as 
a  function  of  only  the  coordinates,  momenta,  and  time: 


T-L  =  'H(qi,pi,t ) 


(11) 
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If  the  Hamiltonian  does  not  contain  time  explicitly,  it  is  a  constant  of  the  motion. 
Taking  the  differentials  of  Equations  10  and  11  yields: 

(YH  =  qdpi  +  pidqi  -  dL(q,  qh  t )  (12) 


and 


on  ,  dn  ,  on  u 

dT-L  =  -7— dqt  +  —  dpt  +  -wrdt 
oqi  dp.i  dt 


respectively.  Equation  12  can  be  simplified  by  expanding  dL  as 


,  dL  dL  dL 
dL  =  —dqi  +  —  dqi  +  —dt 
dqi  dqi  dt 


and  re-writing  in  terms  of  Pi  and  pi  using  Equations  8  and  9 


dL 

dL  =  p^  +  pidqi  +  —  dt 

dt 


(13) 


(14) 


(15) 


which  yields 

dL 

dLi  =  qidpi  -  pidqi  -  —dt  (16) 

dt 

Expressions  for  the  time  derivatives  of  the  coordinates  and  momenta  can  then  be 
written  as  partial  derivatives  of  LL  =  ^(g^p^f)  by  equating  the  first  two  terms  of 
Equations  13  and  16: 

.  dU  .  dn 

Qi  W  5  Pi  W 

uPi  dqi 

Equations  17  are  a  set  of  first-order  differential  equations  of  motion  known  as  Hamil¬ 
ton’s  equations.  A  set  of  g,;  and  Pi  which  satisfy  Hamilton’s  equations  are  said  to  be 
canonical. 

It  is  evident  from  Equation  17  that  if  the  Hamiltonian  can  be  written  indepen¬ 
dent  of  a  coordinate  q,  or  momenta  /y,  then  the  conjugate  momenta  or  coordinate 
will  be  a  constant.  This  demonstrates  that  the  correct  choice  of  coordinates  and 
momenta  can  greatly  simplify  the  dynamics  of  a  system. 
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2. 9  Hamilton- Jacobi  Theory 

Hamilton- Jacobi  Theory  builds  on  the  concepts  of  Hamiltonian  dynamics  in 
that  it  can  be  used  to  find  the  most  convenient  coordinates  and  momenta  for  a 
system  such  that  they  maximize  the  simplicity  of  the  equations  of  motion.  This  is 
done  through  the  use  of  a  generating  function,  F,  which  is  a  function  of  a  mixture 
of  new  and  old  coordinates  and  momenta  and  is  used  to  transform  between  the  new 
coordinates  and  momenta,  Qi:  and  Pt  and  the  old  (qu  pi).  For  a  discussion  on  the 
formulation  of  generating  functions,  see  Wiesel  [35]. 

Hamilton’s  equations  can  be  written  in  the  new  coordinates  and  momenta  as 


Qi 


dlC 

dQi 


(18) 


Where  /C  is  the  system  Hamiltonian  expressed  in  the  new  coordinates  and  momenta 
obtained  from  the  original  Hamiltonian  and  the  generating  function.  The  simplest 
representation  of  a  dynamical  system  would  be  when  Hamilton’s  equations  are  all 
equal  to  0  or,  equivalently,  when  all  of  the  coordinates  and  momenta  are  constants. 
In  this  case,  the  new  Hamiltonian  would  be  exactly  equal  to  0  for  all  time: 

d  F 

1C  =  n(qi,pi,t)  +  —  =  0  (19) 


Where  F  is  the  generating  function.  If  an  F2  generating  function  is  used  (one  written 
in  terms  of  the  old  coordinates,  the  new  constant  momenta  and  time:  F2(qll  Pt,  t)) 
the  old  momenta,  pt  can  be  written  in  terms  of  partial  derivatives  of  F2\ 

8F2 
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and  Equation  19  can  be  re-written  as  a  partial  differential  equation  in  terms  of  the 
generating  function,  here  denoted  S 


dS 


+  —  =  0 


d  q, 


dS 


dt 


(21) 


The  solution  to  Equation  21  gives  the  specific  generating  function  S  which  can  be 
used  to  transform  a  system  to  new  coordinates  and  momenta  which  are  all  constant. 
This  generating  function  is  known  as  Hamilton’s  Principal  Function. 

In  the  special  case  that  the  Hamiltonian  is  independent  of  time,  and  therefore 
a  constant  of  the  motion,  Hamilton’s  Principal  Function  can  be  separated  into  two 
parts:  One  a  function  of  only  the  old  coordinates  and  new  momenta  and  one  a 
function  only  of  time 

5  =  W(qi,Pi)  +  St(t)  (22) 


Where  W  is  known  as  Hamilton’s  characteristic  function.  Using  this  formulation  of 
S,  Equation  21  can  be  written  as 


n 


dSt 

H - - 

dt 


0 


(23) 


which  can  be  separated  to  give 


n 


Pi 


(24) 


If  Hamilton’s  characteristic  function  is  used  as  a  type  2  generating  function, 
the  transformation  equations  between  the  new  and  old  coordinates  can  be  written 
as 

dW  „  dW 

Pi  “a  i  o  r> 

UQi  O  Pi 

where  the  new  momenta,  Ptl  are  all  constant  [35]. 
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The  new  Hamiltonian  can  again  be  written  from  Equation  19,  but  since  the 
generating  function,  W,  is  independent  of  time,  the  new  Hamiltonian  is  equal  to  the 
old  Hamiltonian,  and  since  the  new  momenta  are  all  constant,  the  Hamiltonian  must 
not  include  the  new  coordinates.  That  is 


K.  =  H(Pi) 


(26) 


This  formulation  then  results  in  a  system  in  which  the  new  momenta  are  all  constant 
while  the  new  coordinates  are  all  cyclic.  This  can  be  seen  from  Hamilton’s  equations 


Qi 


dlC 

OR; 


=  v. 


p. 

I  1  1  l 


OK. 

dQi 


(27) 


where  ul  are  constants. 


2.10  Action- Angle  Variables 

In  the  case  above  where  the  new  coordinates  are  cyclic  and  the  motion  of  the 
system  is  periodic,  it  is  possible  to  choose  another  set  of  coordinates  and  momenta 
from  which  the  frequencies  of  the  motion  are  easily  seen.  In  this  case,  it  is  desired  that 
the  coordinates  be  those  quantities  which  are  periodic  and  the  conjugate  momenta  be 
their  associated  momenta.  A  set  of  coordinates  and  momenta  of  this  type  are  known 
as  Action- Angle  variables.  The  momenta,  or  the  action  variables  J,  can  be  calculated 
from  the  old  coordinates  and  momenta  using  the  principle  of  least  action  [11], 


J  = 


(28) 


Where  the  integration  is  carried  out  over  a  complete  period  of  the  motion.  The 
conjugate  coordinates,  or  angles  w,  can  then  be  calculated  by  first  writing  Hamilton’s 
characteristic  function 

W  =  W{q ,  J)  (29) 
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and  using  Equation  25  to  find 


(30) 


dW 
W  ~  <).l 

The  choice  of  w  and  J  as  system  coordinates  and  momenta  is  beneficial  because  it 
can  be  shown  that  the  time  rates  of  change  of  the  coordinates,  w,  are  exactly  the 
frequencies  associated  with  the  periodic  motion  of  the  system  [11], 

2.11  KAM  Theorem 

Kolmogorov,  Arnol’d,  and  Moser  Theorem,  or  KAM  Theorem,  was  originally 
developed  in  the  1950s  as  an  attempt  to  solve  a  type  of  problem  first  come  upon 
in  the  field  of  celestial  mechanics  in  the  1700s  [31].  This  problem,  now  known  as 
the  three-body  problem  (3BP),  was  posed  by  Newton  who  had  written  differential 
equations  describing  the  interaction  of  multiple  massive  bodies  under  the  influence 
of  gravity.  This  problem  was  shown  to  have  a  closed-form  solution  if  there  were  only 
two  bodies  (the  two-body  problem,  (2BP))  but  if  a  third  body  was  introduced,  no 
closed  form  solution  could  be  found.  If  two  of  the  bodies  are  much  less  massive  than 
the  third,  the  system  can  be  expressed  as  the  integrable  2BP  with  a  perturbation  due 
to  the  third  body.  This  problem  was  then  classically  ‘solved’  using  series  expansions, 
but  convergence  was  not  possible  due  to  the  appearance  of  small  divisors  caused  by 
resonances  in  the  Hamiltonian  dynamics.  In  1954,  A.N.  Kolmogorov  suggested  two 
ideas  which  are  central  to  the  KAM  technique  [31]; 

1.  Linearize  the  problem  about  an  approximate  solution  and  solve  the  linearized 
problem 

2.  Improve  the  approximate  solution  by  using  the  linearized  problem  solution  as 
the  basis  of  a  Newton- Rhapson  method  argument 

These  central  ideas  were  built  upon  by  V.  Arnol’d  [1]  and  J.  Moser  in  the  1950s  and 
1960s  and  came  to  be  known  as  KAM  theorem. 
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In  words,  KAM  theorem  states  that  for  a  integrable  system  in  which  the  mo¬ 
menta  and  forces  are  invariant  -  that  is,  a  Hamiltonian  system  -  subject  to  small 
smooth  perturbations  from  conservative  forces,  many  of  the  solutions  for  the  unper¬ 
turbed  system  are  also  solutions  to  the  perturbed  system  with  small  changes  [35]. 
The  perturbed  system  contains  one  or  more  action-angle  variables  such  that  the 
Hamiltonian  can  be  written  as 

H(J,w)  =  h(J)  +  ef(J,w)  (31) 

where  J  and  w  are  the  action-angle  variables,  h  is  the  unperturbed  Hamiltonian, 
/  is  the  perturbing  function  and  e  is  the  small  perturbing  parameter.  Therefore, 
for  e  =  0  in  Equation  31,  the  Hamiltonian  reduces  to  the  initial,  integrable  system. 
Solutions  to  the  system  then  can  be  shown  to  have  the  characteristic  of  returning  to 
their  initial  position  if  one  angular  coordinate  is  incremented  by  an  integer  multiple 
of  some  characteristic  angular  period  while  the  others  are  held  fixed.  This  type  of 
solution  is  defined  as  a  torus.  The  full  mathematical  description  and  proof  of  KAM 
theorem  can  be  found  in  [25]  for  example. 

A  torus  can  be  thought  of  geometrically  as  the  product  of  N  circles,  where  N 
is  then  the  dimension  of  the  torus.  For  example,  a  one-dimensional  torus  is  simply 
a  circle,  while  a  two-dimensional  torus  is  the  product  of  two  circles;  Think  of  the 
surface  formed  by  revolving  a  smaller  circle  around  the  perimeter  of  a  larger  circle,  or 
a  donut.  An  N  dimensional  torus  exists  in  2N  dimensional  phase  space.  Continuing 
with  the  product  of  circles  analogy,  the  2N  dimensions  would  be  the  radius  of  each 
circle,  the  actions,  J,  and  the  angles,  w,  would  be  the  angles  between  vectors  from 
the  center  of  each  circle  to  the  desired  location  on  the  circle  and  some  reference  line. 
Note  that  the  actions  are  all  constant.  Figures  6  and  7  show  pictures  of  one-  and 
two-dimensional  tori  and  their  action-angle  variables. 
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Figure  6:  Diagram  of  a  One-Dimensional  Torus 


Top  View  Cut-Away  View 

Figure  7:  Diagram  of  a  Two-Dimensional  Torus 


Higher  dimensional  tori  are  hard  to  visualize  as  they  cannot  be  explicitly  drawn 
in  three  physical  dimensions.  However,  an  example  of  a  three-dimensional  torus  can 
be  thought  of  in  the  following  way.  Consider  a  satellite  traveling  around  an  orbit 
with  both  a  precessing  argument  of  perigee  and  a  regressing  node.  If  both  u  and  C2 
are  held  constant,  and  the  mean  anomaly  is  advanced  by  n2n  where  n  is  any  integer, 
the  satellite  will  return  to  the  exact  position  at  which  it  started.  In  fact,  if  any  two 
of  the  three  coordinates,  (cu,  C2,  M)  are  held  fixed  while  the  third  is  incremented  by 
an  integer  multiple  of  2n,  the  satellite’s  position  will  remain  unchanged.  Therefore, 


22 


the  satellite  can  be  said  to  be  traversing  a  3-dimensional  torus  in  6-dimensional 
phase  space,  with  the  dimensions  being  the  three  angular  coordinates,  and  the  three 
conjugate  action  momenta  obtained  from  the  system  Hamiltonian  and  Hamilton’s 
equations. 

2.12  Reference  Frames 

The  discussions  of  Sections  2.8  through  2.11  can  be  applied  to  the  equations  of 
motion  for  an  Earth-orbiting  satellite.  There  are  several  coordinate  frames  that  can 
be  useful  in  describing  the  motion  of  an  Earth-orbiting  satellite.  Three  frames  that 
will  be  used  in  this  work  are  the  Perifocal  (PF)  Frame,  the  Earth-Centered  Inertial 
(ECI)  frame,  and  the  Earth-Centered-Earth-Fixed  (ECEF)  frame.  These  frames  are 
defined  in  Table  2. 


Table  2:  Coordinate  Frame  Definitions 


Coordinate 

frame 

Type 

Origin 

1-axis 

3-axis 

2-axis 

PF 

Inertial 

Earth  CoM 

Toward  perigee 

Along  angular 
momentum  vector 

Completes  right- 
handed  frame 

ECI 

Inertial 

Earth  CoM 

Toward  vernal 
equinox 

Earth’s  axis 
of  rotation 

Completes  right- 
handed  frame 

ECEF 

Non- 

inertial 

Earth  CoM 

Toward  prime 
meridian  in 
equatorial  plane 

Earth’s  axis 
of  rotation 

Completes  right- 
handed  frame 

2.13  Earth  Satellite  Dynamics 

Let  the  position  of  a  satellite  be  denoted  in  rectangular  coordinates  in  the 
ECEF  frame  as  x,  y,  z.  The  inertial  velocity,  resolved  along  the  ECEF  axes  can  then 
be  written  as 

v  =  y  +  UJ&X  (32) 

i 
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Where  cu©  is  the  rotational  rate  of  the  Earth  [33].  The  kinetic  energy,  per  unit  mass 
of  the  satellite,  can  be  written  as 

T  =  -v  ■  v  —  -  [(i;  —  t o^yY  +  (y  +  u>®x)~  +  z2}  (33) 

and  the  canonical  momenta,  again  per  unit  mass,  are  just  the  velocity  components 
from  Equation  32.  The  potential  energy  for  the  system  can  be  written  per  unit 
mass  as  the  Earth’s  geopoential  given  in  Equation  5.  The  Hamiltonian  can  then  be 
written  from  Equation  10,  with  the  (/,;  being  written  in  terms  of  the  momenta  from 
Equation  32. 

ft  =  \  (pI  +  P2y  +  Pi)  +  (VPx  ~  XPy)  +  V  (34) 

where  V  is  the  Geopotential  from  Equation  5.  Pi  is  independent  of  time  and  is 
therefore  a  constant  of  the  motion.  Note  that  the  coordinates  in  Equation  34  are 
rectangular  while  the  Geopotential  was  expressed  in  spherical  coordinates.  The 
transformation  between  the  two  sets  of  coordinates  is  given  by 


r  = 
sin  5  = 


\J x2  +  y2  +  z 
z 

\/x2  +  y2 


2 


tan  A 


y 

X 


(35) 


where  r  is  the  scalar  radius,  5  is  the  latitude,  and  A  is  the  east  longitude. 

The  motion  of  an  Earth-orbiting  satellite  is  obviously  periodic.  Wiesel  has 
shown  how  the  basis  frequency  set  of  an  earth-orbiting  satellite  can  be  described  by 
analyzing  its  position  as  a  function  of  time  as  follows  [33].  In  the  Perifocal  reference 
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frame  the  satellite  position  can  be  written  as 

r cos(u ) 

rsin(z/)  (36) 

0 

where  r  is  the  scalar  orbital  radius,  and  v  is  the  true  anomaly,  v  can  be  written  in 
terms  of  the  mean  anomaly,  M,  and  therefore  position  can  be  written  as  a  periodic 
function  of  the  mean  anomaly,  rpf  =  rpf(M).  This  position  can  be  transformed  to 
the  ECEF  frame  using  simple  1  and  3-axis  rotation  matrices 

reCef  =  R3(-9)R3(-tt)R1(-i)R3(-to}rpf(M)  (37) 

Where  9  is  the  Greenwich  sidereal  time,  is  the  right  ascension  of  the  ascending 
node  (RAAN),  i  is  the  inclination,  and  to  is  the  argument  of  perigee.  9 ,  and  u 
can  be  written  as  functions  of  time  as: 

9  =  9q  +  c U($t 

Q  =  Q0  +  Clt  (38) 

UJ  =  LOo  +  LOt 

where  the  0  subscript  denotes  value  at  epoch.  Plugging  these  values  into  Equation  37 
and  re-arranging  yields 

fecef  =  R3(90  -  fi0)-R3((we  -  n)t)i?i(-i)i?3(-w0)i?3(-d;t)rp/(M)  (39) 

From  this  form,  it  is  apparent  that  the  position  vector  in  the  ECEF  frame  will  be 

dependent  on  three  periodic  terms  which  become  the  basis  frequency  set  for  a  KAM 
torus:  M,  0%  —  f2,  and  —  Co.  These  terms  can  be  written  in  terms  of  the  other  COEs 
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and  the  earth’s  gravitational  field,  through  J2,  as 


Cm  =  M  =  yf$ 
0J2 =  n  —  cu0  — 


3J2-R0 


2a2(l— e2)1 


(hi 


sm2fi) 


1) 


3^JU2R2 


2a3-5(l— e2)2  COSOV  W© 


cn3  =  cfi  = 


2a3-5(l— e2)2 


(|  sin2(i)  -  l) 


(40) 


Where  e  is  the  orbital  eccentricity,  a  is  the  semi-major  axis,  /i  is  the  gravitational 
parameter,  and  Pe  is  the  radius  of  the  earth.  For  a  derivation  of  M,  D,  and  cfi, 
see  Danby  [8].  Note  the  cm  and  cm  have  switched  sign  from  Equation  39  to  Equa¬ 
tion  40.  This  is  done  to  keep  the  torus  frequencies  consistent  with  normal  Keplarian 
mechanics;  That  is,  the  line  of  apsides  precesses  while  the  line  of  nodes  regresses. 

These  frequencies  can  then  be  used  to  form  the  action-angle  variables  for  the 
system.  Specifically,  the  angles  are  the  linear  coordinates  with  time  derivatives  equal 
to  the  basis  frequencies: 

Qi  =  M 

Q2  =  cnet  (41) 

Q3  —  ^ 

The  conjugate  action  momenta  can  be  calculated  from  Equation  28  and  have  been 
shown  to  be  approximately  equal  to  the  Delaunay  momenta  [34], 


Pi  =  y/JW, 

P2  =  y[ga\J\  —  e2  cos i  (42) 

P3  =  yfgaVl  -  e2 


The  Hamiltonian,  Equation  34,  can  then  be  re-written  as  a  function  of  only  the  new 
momenta  as 


K.  = 


/P 


—  CUm  + 


^J2R%  (Pi  -  3 Pi) 


(43) 


2P2  w  4  PfPf 

[34],  The  action-angle  variables  in  Equations  41  and  42  can  then  be  used  to  express 
the  motion  of  an  Earth-orbiting  satellite  as  a  KAM  torus. 
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2.14  Application  to  Celestial  Mechanics 

As  noted  in  Section  2.11,  KAM  theorem  can  be  applied  to  lightly  perturbed 
Hamiltonian  systems.  Since  the  initial  theorem  was  developed,  some  effort  has  been 
made  to  apply  the  theorem  to  celestial  mechanics. 

Celletti  has  done  work  to  apply  KAM  theorem  to  a  number  of  celestial  mechan¬ 
ics  problems  including  the  restricted  3BP,  the  planetary  N-Body  problem,  and  the 
Spin-Orbit  problem  in  which  a  rigid,  tri-axial  satellite  orbits  a  central  body  and  the 
orbital  revolution  and  rotational  motion  of  the  satellite  are  coupled  [4] .  In  analyzing 
the  specific  restricted,  circular,  planar  3BP  of  the  asteroid  12  Victoria  orbiting  the 
Sun  and  perturbed  by  Jupiter’s  gravitational  held,  Celletti  and  Chiercha  demon¬ 
strated  an  application  of  an  iso-energetic  KAM  method  in  which  invariant  ‘trapping’ 
tori  are  constructed  on  an  energy  level  similar  to  the  osculating  Keplarian  motion 
of  the  asteroid.  These  trapping  tori  then  bound  the  values  of  the  action  variables  of 
the  asteroid  [5]. 

Robutel  demonstrated  an  application  of  KAM  theorem  to  the  planetary  3BP 
using  canonical  heliocentric  variables  [26].  This  choice  of  variables  simplifies  the 
system  Hamiltonian  and  allows  the  perturbing  function  to  be  written  in  a  compact 
form.  The  problem  is  reduced  to  a  four  degree  of  freedom  system,  resulting  in  a 
KAM  torus  in  4-dimensional  phase-space.  This  torus  is  shown  to  be  stable  for  all 
time  for  small  planetary  masses  in  orbits  with  small  eccentricity  and  inclination. 

McGill  and  Binney  developed  method  of  generating  tori  for  a  general  gravita¬ 
tional  system  using  a  type-2  generating  function  to  map  between  the  action-angle  co¬ 
ordinates  for  a  well-known  potential  and  those  of  the  system  being  investigated  [22] . 
Their  method  utilizes  a  non-linear  least-squares  technique  to  determine  coefficients 
for  the  generating  function  and  distorts  a  ‘toy  torus’  from  the  well-known  system 
into  a  torus  for  the  system  of  interest.  They  applied  their  technique  to  find  a  torus 
for  Keplarian  motion  starting  from  the  torus  of  a  degenerate  harmonic  oscillator. 
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2.15  Application  to  Earth- Orbiting  Satellites 

To  date,  the  only  work  found  related  to  the  application  of  KAM  theorem  to 
Earth-orbiting  satellites  has  been  done  by  Wiesel  and  his  Masters  and  PhD  students. 

Wiesel’s  first  article  on  the  subject  showed  that  the  KAM  theory  could  best 
be  applied  to  earth  orbiting  satellites  when  their  motion  was  expressed  in  the  Earth- 
Centered-Earth-Fixed  (ECEF)  coordinate  frame  [33].  In  this  frame,  the  Earth’s 
gravitational  held  is  nearly  constant  and  the  the  system  momenta  can  be  written 
simply  as  the  time  derivatives  of  the  coordinates,  that  is,  as  the  velocity  components 
as  demonstrated  in  Section  2.13.  The  Hamiltonian  can  then  be  written  indepen¬ 
dent  of  time  in  terms  of  only  the  momenta  and  the  gravity  held  and  is  therefore 
an  integral  of  the  motion.  Using  orbital  position  data  obtained  from  a  one-year 
numerical  integration,  Wiesel  extracted  the  frequency  content  of  the  orbit  using  a 
Fast  Fourier  Transform  (FFT)  and  showed  that  for  nearly  circular  orbits  with  mod¬ 
erate  inclination,  the  frequency  spectrum  exhibited  clearly  defined  peaks  which  were 
linear  combinations  of  only  three  discrete  basis  frequencies,  showing  that  their  mo¬ 
tion  might  lie  on  a  KAM  torus.  For  orbits  with  inclinations  near  90  degrees,  the 
frequency  spectra  were  not  as  clean,  exhibiting  nearly  chaotic  behavior.  In  order  to 
more  precisely  define  the  basis  frequencies,  a  method  was  used  based  on  the  work 
of  Laskar  [19].  Having  obtained  the  basis  frequency  set,  the  orbital  data  was  fit 
to  a  Fourier  series  in  the  basis  frequencies.  This  Fourier  series  contains  the  torus 
coordinates  and  is  used  to  transform  from  torus  coordinates  to  physical  cartesian 
coordinates.  A  version  of  this  process  was  used  for  the  current  effort  and  will  be 
described  in  Section  3.5. 

In  subsequent  work,  Wiesel  further  analyzed  the  frequency  spectrum  of  a  earth¬ 
orbiting  satellite  from  orbital  data  obtained  using  a  simplified  version  of  SGP4  [34], 
This  analysis  showed  that  the  most  prominent  spectral  lines  occur  in  clusters  around 
multiples  of  the  anomalistic  frequency,  M,  separated  by  combinations  of  the  rota¬ 
tional  rate  of  the  Earth,  the  nodal  regression  rate  and  the  rate  of  perigee  precession. 


Wiesel  then  compared  the  trajectories  obtained  from  a  KAM  torus,  constructed  by 
the  method  described  by  in  [33],  with  tori  obtained  using  the  simplified  SGP4  model 
and  the  2BP.  The  trajectories  from  these  tori  were  then  compared  to  a  numerically 
integrated  trajectory.  The  KAM  torus  was  accurate  to  within  30  meters  over  six 
months  while  the  SGP4  torus  was  accurate  only  to  within  40  kilometers. 

Derbis  analyzed  precise  orbital  data  from  26  GPS  satellites  using  a  Fast  Fourier 
Transform  (FFT)  to  extract  approximate  basis  frequencies  [9].  These  frequencies 
were  compared  to  those  obtained  from  numerically  integrated  orbital  data  and  the 
Laskar  frequency  method.  The  two  methods  yielded  similar  but  not  exact  duplicate 
frequency  sets.  Her  research  also  illustrated  two  difficulties  in  applying  KAM  theo¬ 
rem  to  artificial  satellites.  If  pairs  of  basis  frequencies  are  near-integer  multiples  of 
each  other  (nearly  commensurate),  the  basis  set  becomes  hard  to  determine  because 
the  spectral  lines  lie  nearly  on  top  of  each  other.  In  Derbis’  case,  this  difficulty 
presented  itself  in  oj2  ~  2ui  because  the  GPS  constellation  has  an  orbital  period  of 
12  hours,  or  one-half  the  rotational  period  of  Earth.  The  other  difficulty  was  seen 
when  analyzing  data  from  the  oldest  GPS  satellite  which  uses  only  3  reaction  wheels 
for  attitude  control  and  must  periodically  do  momentum  dumps.  These  momentum 
dumps  result  in  nearly-impulsive  A v's  which,  although  small,  do  not  fit  the  KAM 
theorem  criteria  of  small  smooth  perturbations.  As  a  result  of  this,  the  spectral  map 
of  this  satellite  exhibited  significant  noise  and  revealed  basis  frequencies  which  dif¬ 
fered  from  the  other  GPS  satellites.  This  shows  that  KAM  theorem  must  be  applied 
carefully  to  operational  satellites  that  maneuver  from  time-to-time  as  a  large  enough 
maneuver  will  move  the  satellite  off  of  the  torus.  Current  research  is  being  done  on 
the  behavior  of  motion  near  a  KAM  torus  in  the  hope  of  determining  a  way  to 
account  for  this. 

Little  analyzed  orbital  data  from  the  Gravity  Recovery  and  Climate  Experi¬ 
ence  (GRACE)  and  Jason-1  satellites  and  determined  their  basis  frequencies  using 
the  Laskar  method  [20].  The  predicted  position  from  the  resulting  torus  was  then 
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compared  to  the  actual  position  data  to  determine  its  accuracy.  Based  on  the  growth 
rate  of  the  position  error,  the  basis  frequencies  were  adjusted  which  resulted  in  a 
torus  which  provided  significant  improvement  in  position  accuracy.  Using  this  pro¬ 
cess,  KAM  torus  position  data  for  the  Jason-1  satellite  was  shown  to  be  accurate  to 
about  1  kilometer  over  15  days,  while  the  torus  for  the  GRACE  satellite  was  accu¬ 
rate  to  only  approximately  15  kilometers  over  the  same  time  period.  The  residual 
growth  in  the  GRACE  data  was  quadratic  which  was  taken  to  be  a  result  of  air  drag. 
Because  air-drag  is  a  non- conservative  force,  it  cannot  be  readily  incorporated  into 
the  KAM  theorem.  This  is  also  a  topic  of  current  research. 

Craft  performed  an  analysis  of  the  effect  of  the  number  of  terms  in  the  torus 
Fourier  series  to  position  accuracy,  measured  against  a  numerically  integrated  or¬ 
bit  [7].  He  showed  that  the  ideal  number  of  series  terms,  with  respect  to  accuracy 
and  computational  burden,  is  between  1500  and  2000.  For  his  test  cases,  this  number 
of  terms  resulted  in  an  rms  error  in  the  position  components  on  the  order  of  10_2km. 
Craft  also  studied  the  applicability  of  KAM  theorem  to  satellite  formation  flight  [7] . 
In  this  analysis,  a  torus  was  constructed  for  a  chief  satellite  and  a  deputy  satellite 
was  placed  on  the  same  torus  at  an  offset  of  one  or  more  of  the  torus  angle  coor¬ 
dinates.  This  resulted  in  a  position  difference  between  the  chief  and  deputy  which 
was  oscillatory  in  time  and  contained  a  small  secular  drift  rate  proportional  to  the 
initial  displacement. 

Bordner  attempted  to  fit  KAM  tori  to  high-precision  orbital  data  for  GPS 
satellites  using  a  variety  of  spectral  methods  based  on  those  developed  by  Laskar  and 
Wiesel,  but  was  unable  to  produce  a  torus  with  suitable  accuracy  [2],  He  encountered 
the  same  difficulties  with  near-commensurate  frequencies  as  Derbis,  and  also  found 
cu3  very  hard  to  identify  using  spectral  methods  because  its  period  was  very  long 
compared  to  the  timespan  of  the  data  being  analyzed.  Another  method  of  developing 
the  torus  Fourier  series  was  attempted,  this  time  using  a  least  squares  approach  to 
fit  the  coefficients  of  the  series  after  obtaining  the  basis  frequencies.  This  method 
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yielded  improved  accuracy,  but  still  had  errors  in  excess  of  acceptable  levels  for  the 
GPS  constellation. 

Bordner  also  developed  methods  for  torus  construction  from  numerically  in¬ 
tegrated  trajectories.  The  most  successful  method  was  shown  to  be  one  in  which 
clusters  of  peaks  in  the  frequency  power  spectrum  were  decomposed  simultaneously. 
This  method  allowed  CU3  to  be  identified  more  easily  as  the  separation  between  the 
dominant  peak  and  its  flanking  peaks  in  each  cluster.  Using  this  technique,  position 
errors  were  shown  to  drop  to  10s  of  meters  over  a  period  of  6  months  for  certain 
low-earth  orbits  [2], 

2.16  Contribution  of  Current  Work 

The  current  work  will  build  on  the  results  described  above  in  Section  2.15 
in  a  number  of  ways.  First,  an  attempt  will  be  made  to  identify  KAM  torus  basis 
frequencies  from  purely  observational  data  in  the  form  of  TLEs.  This  is  in  contrast  to 
the  current  method  of  obtaining  the  basis  frequencies  from  a  numerically  integrated 
trajectory.  A  Matlab  script  will  be  developed  to  read  in  and  format  TLE  data 
to  enable  basis  frequency  identification  through  least-squares  curve  fitting.  The 
script  will  also  read  in  and  format  additional  TLE  data,  such  as  the  bstar  drag 
term,  which  could  be  used  to  build  on  the  current  effort  in  future  work.  This  work 
will  also  demonstrate  a  new  technique  of  making  small  changes  to  a  torus’  basis 
frequencies  by  first  calculating  a  torus  momenta  offset  corresponding  to  a  known 
frequency  offset  and  then  translating  that  momenta  offset  to  a  change  in  initial 
velocity.  This  frequency  matching  process  can  be  utilized  in  an  iteration  with  the 
current  torus  construction  algorithm  to  obtain  a  torus  with  basis  frequencies  nearly 
exactly  equal  to  the  desired  values.  Finally,  this  work  includes  a  detailed  analysis  of 
the  frequency  content  of  torus  position  error  data  which  can  be  used  going  forward 
to  refine  the  torus  Fourier  series  construction  algorithms  to  produce  tori  with  more 
accurate  position  prediction  capability. 
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III.  Methodology 

3.1  Test  Object  Selection  and  Data  Gathering 

The  TLEs  used  for  this  effort  were  obtained  from  the  website  http: / / celestrak.com. 
This  site  allows  the  user  to  request  historical  TLEs  for  any  object  in  the  satellite  cat¬ 
alog.  As  discussed  in  Section  2.15,  previous  research  has  demonstrated  that  certain 
types  of  orbits  cause  difficulties  in  applying  KAM  theorem  to  Earth  satellites.  These 
difficulties  were  taken  into  account  in  selecting  test  objects  for  this  effort.  Table  3 
lists  each  known  issue  and  describes  how  each  was  accounted  for  in  selecting  the 
current  test  objects. 


Table  3:  Test  Object  Selection  Criteria 


Issue 

Object  Selection  Criteria 

Air  Drag 

Choose  objects  orbiting  above  300km  altitude  that 
have  been  in  stable  orbits  for  a  long  period  of  time. 

Station  Keeping 

Choose  either  non-operational  objects  or  operational 
objects  which  use  only  reaction  wheels  for  attitude  control. 

Near- Commensurate  Frequencies 

Choose  objects  such  that  the  Earth’s  rotational  period  is 
not  nearly  an  integer  multiple  of  the  orbital  period. 

Near-Polar  Inclination 

Choose  objects  with  inclination  below  critical. 

This  analysis  led  to  the  choices  of  three  spent  rocket  bodies  and  the  HST  as 
test  subjects  for  this  effort.  General  information  on  the  approximate  orbits  of  each 
of  these  satellites  is  in  Table  4. 


Table  4:  Orbital  Information  for  Test  Satellites 


Name 

Catalog  Num. 

Launch  Date 

Period 

[min] 

Inclination 

[deg] 

Apogee  altitude 
[km] 

Perigee  altitude 
[km] 

Eccentricity 

Hubble  Space  Telescope 

20580 

4/24/1990 

95.93 

28.47 

566 

561 

3.836E-4 

Thor  Ablestar  Rocket  Body 

59 

10/4/1960 

106.44 

28.25 

1203 

921 

1.893E-2 

Delta  1  Rocket  Body  No.l 

341 

7/10/1962 

157.52 

44.77 

5619 

949 

2.414E-1 

Delta  1  Rocket  Body  No. 2 

8133 

8/27/1975 

95.21 

25.30 

700 

357 

2.753E-2 

A  Matlab  script  was  written  to  read  a  text  hie  containing  multiple  TLEs.  This 
raw  data  was  then  converted  to  units  and  format  useful  for  further  analysis.  All 
angular  data  in  TLEs  is  reported  in  degrees,  from  0  to  360.  These  values  were 
converted  to  radians,  and  the  2ir  jumps  were  eliminated,  resulting  in  smooth  plots 
of  RAAN  and  argument  of  perigee  as  functions  of  time.  The  mean  anomaly  data 
was  merged  with  the  revolution  number  data  to  form  a  continuous  plot  of  M  vs. 
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time.  Time  in  TLEs  is  represented  as  a  number  between  0  and  367  which  gives  the 
day  (and  partial  day)  of  the  year,  starting  from  0000  UT  31  Dec.  For  example,  a 
epoch  time  of  09001.000000000  corresponds  to  0000  UT  on  01  Jan  2009,  while  an 
epoch  time  of  09000.00000000  corresponds  to  an  epoch  time  of  0000  UT  on  31  Dec 
2008.  All  times  are  measured  in  mean  solar  days  (24  hour  days)  rather  than  sidereal 
days.  The  Matlab  script  converted  this  raw  data  to  a  continuous  timescale,  starting 
at  zero,  and  changed  units  from  mean  solar  days,  to  sidereal  days,  to  canonical  time 
units  (TUs)  as  was  necessary. 


3.2  TIE  Frequency  Identification 

After  the  data  was  read-in  and  formatted,  further  analysis  was  done  to  extract 
the  characteristic  frequencies  of  the  orbit,  that  is: 


dM  dil  du 
~dT  '  ~dt  '  ~dt 


(44) 


To  do  this,  a  second  order  curve-fit  was  accomplished  using  a  least  squares 
technique  as  follows: 


Let  the  data  to  be  fit  be  represented  as  a  vector  d  while  the  time  corresponding 
to  each  data  point  is  grouped  into  a  vector  t.  Then,  the  second  order  curve  fit  will 
be  of  the  form 


d  = 


T  (if  +  —02^ 


+2 


(45) 


where  at  are  the  curve-fit  coefficients.  Next,  define  a  matrix  T  as 


T  = 


dd 

dai 


1  t 


If1 

2L 


(46) 
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where  1  is  a  column  vector  of  the  same  length  of  t  containing  all  ones.  Now,  the 
curve-fit  coefficients  can  be  solved  for 

a  =  (TtQ_1T)  _1  TTQ~ld  (47) 

where  Q  is  assumed  to  be  the  identity  matrix  for  now  since  the  individual  data  point 
covariances  are  not  known.  The  curve  fit  is  then  given  by 

/  =  Td  (48) 

and  residuals,  r,  can  be  calculated: 

r  =  d  —  f  (49) 

To  determine  the  accuracy  of  the  curve-fit,  the  covariance  matrix,  P,  was  calculated 

P=(TTQ-1Ty1f0  (50) 

where  f o  is  the  average  squared  residual,  which  we  take  to  be  approximately  equal 
to  the  standard  deviation  squared,  a2 

1  N 

=  ~  a2  (51) 
i= 1 

and  N  is  the  number  of  data  points. 

For  the  HST  for  example,  these  curve  fits  came  out  to: 

M  =  1.166155  x  10 ~5t2  +  9.425866  x  101!  +  3.487986 
n  =  -3.228504  x  10"8t2  -  1.137980  x  10 ~H  +  7.145640  x  10"1  (52) 

u  =  -4.202363  x  10 ~H2  +  1.853485  x  10_1t  +  2.714583 
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With  time  in  units  of  mean  solar  days.  Note  that  the  coefficients  have  been  rounded 
and  more  significant  digits  were  carried  in  the  actual  calculations.  As  expected,  the 
second  order  terms  in  the  curve-fits  are  much  smaller  than  the  lower  order  terms. 
These  second  order  terms  appear  due  to  the  presence  of  air-drag  and  other  low-order 
effects.  In  forming  the  initial  KAM  tori,  the  desired  basis  frequencies  were  based  on 
the  first  order  curve-fit  coefficients  and  second  order  effects  were  accounted  for  later. 

3.3  Canonical  Units 

In  many  of  the  calculations  performed,  canonical  units  of  distance  units  (DUs) 
and  time  units  (TUs)  were  used.  The  values  used  for  these  quantities  were 

1  DIJ  =  6378.135  km 

(53) 

1  TU  =  13.44686457  min 

3.4  Position  and  Velocity 

In  order  to  numerically  integrate  the  trajectory  and  obtain  a  data  set  from 
which  to  build  a  torus,  an  initial  position  and  velocity  is  necessary.  This  was  ob¬ 
tained  from  the  TLEs  using  the  SGP4  algorithms  developed  by  Vallado  in  C++  and 
subsequently  translated  to  Matlab.  This  code  is  available  for  download  on  the  in¬ 
ternet  [29].  The  SGP4  algorithm  outputs  position  and  velocity  in  the  True  Equator 
Mean  Equinox  (TEME)  reference  frame,  which  is  nearly  the  ECI  frame.  As  discussed 
in  Section  2.15,  position  and  velocity  data  are  needed  in  the  ECEF  frame.  This  was 
accomplished  using  a  modified  version  of  the  Matlab  routine  teme2ecef,  also  from 
Vallado,  to  output  position  and  velocity  in  the  Pseudo  Earth-Fixed  (PEF)  frame 
which  is  equivalent  to  the  ECEF  frame  described  in  Section  2.12.  As  an  example, 
Figure  8  shows  the  position  of  the  Thor  Rocket  Body  at  the  epoch  time  of  each  TLE. 
While  this  data  does  not  appear  to  look  like  an  orbit,  it  is  important  to  recognize 
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that  the  data  points  are,  on  average,  about  12  hours  apart  and  therefore  do  not  show 
the  continuous  trajectory. 


SGP4  Position 


xio5 


X 


Figure  8:  SGP4  Position  at  each  TLE  Epoch  for  Thor  Rocket  Body 


3.5  Torus  Construction 

3.5.1  Numerical  Integration.  The  initial  position  and  velocity  vectors 
obtained  from  the  first  TLE  using  SGP4  were  input  into  a  numerical  integration 
routine  developed  by  Wiesel  which  calculates  a  trajectory  taking  into  account  the 
first  20  terms  in  the  Earth’s  geopotential  from  EGM96  described  in  Section  2.7. 
The  routine  is  a  4th  order  predictor-corrector  algorithm  which  runs  first  backward 
and  then  forward  in  time.  Error  checking  is  done  by  calculating  the  Hamiltonian 
given  in  Equation  34  at  each  time  step  and  ensuring  its  value  does  not  change.  For 
the  current  analysis,  the  orbits  were  integrated  forward  and  backward  for  6  months, 
resulting  in  one  year’s  worth  of  position  data.  For  all  cases,  error  in  the  Hamiltonian 
did  not  exceed  CHO-13. 
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3.5.2  Torus  Frequency  Identification.  The  numerical  data  was  then  an¬ 
alyzed  to  identify  the  orbit’s  fundamental  frequencies  given  in  Equation  40.  To 
accomplish  this,  a  modified  Laskar  frequency  algorithm  developed  by  Wiesel  was 
used  [35].  A  finite  Fourier  transform  of  the  form 

^  J  Q(t)etxp(t/T)dt  (54) 

was  performed  on  the  numerically  integrated  data  where  to  is  frequency,  q  is  the 
physical  coordinate  (x,  y,  or  z),  t  is  time,  T  is  the  total  integration  time,  and  Xv  is 
a  Hann  window  function  given  by 

<55) 

where  p  is  a  parameter  which,  when  increased,  helps  to  widen  the  central  peak  in 
4>(oo)  while  reducing  the  magnitude  of  any  side  lobes.  Care  must  be  given  when 
increasing  p  however,  for  if  peaks  in  <f{to)  are  close  together,  a  value  of  p  which  is 
too  high  can  cause  the  larger  peak  to  obscure  the  smaller  one.  For  this  work  a  value 
of  p  =  2  was  used. 

To  identify  the  basis  frequency  set,  the  power  spectrum,  P  =  |0|2,  was  com¬ 
puted.  Previous  research  by  Wiesel  has  shown  that  the  power  spectra  for  Earth¬ 
orbiting  satellites  contain  multiple  clusters  of  three  peaks  centered  around  integer 
multiples  of  the  first  basis  frequency  (M  or  oof)  [34],  Figure  9  shows  the  power  spec¬ 
trum  from  the  Thor  Rocket  Body  and  illustrates  the  first  two  of  these  peak  clusters, 
which  are  the  largest  in  magnitude  and  are  centered  around  uo\  and  2tj\ .  The  sepa¬ 
ration  between  the  peaks  in  each  cluster  is  the  second  basis  frequency,  u>2.  Note  that 
the  major  peaks  in  the  X  and  Y  spectrum  are  identical  while  the  Z  peaks  are  offset. 
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Partial  Frequency  Power  Spectrum  for  Thor  Rocket  Body 


Figure  9:  Thor  Rocket  Body  Frequency  Power  Spectrum  Illustrating  Peak  Clusters 
Around  uq  and  2u\ 


The  values  of  each  of  these  peaks  is,  as  linear  combinations  of  uq,  u2,  and  oq: 


xla  =  [1  1  If 
*1  =  [1  0  if 

Xib  =  [1  -  1  if 

*2a  =  [2  1  If 
^2  =  [2  0  if 

x2b  =  [2  -  1  if 


(56) 


Local  maxima  in  the  power  spectrum  are  then  searched  for  at  these  locations  using  a 
Newton- Rhapson  method,  starting  from  the  approximate  peak  locations  calculated 
from  Equation  56  and  using  the  J2  frequencies  given  in  Equation  40.  Having  found 
the  actual  locations  of  the  three  peaks  with  the  largest  amplitude,  a  simple  linear 
system  can  be  solved  for  the  actual  values  of  the  basis  frequencies:  cui,  u 2,  and  uq. 
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3.5.3  Fourier  Series  Construction.  Having  obtained  the  torus  basis  fre¬ 
quency  set,  in  order  to  describe  the  satellite  motion  as  a  KAM  torus,  the  physical 
coordinates,  q%  must  be  expressed  as  a  Fourier  Series  in  the  basis  frequency  set  and 
torus  coordinates.  This  Fourier  series  can  be  written  as 

q  =  Yl  cos  3 -  Q  +  Sj  sin  3  '  Q)  (57) 

1 

where  j  is  a  vector  summation  index,  each  vector  entry  identifying  a  spectral  peak 
as,  for  example,  in  Equation  56,  C  and  S  are  the  yet  to  be  solved  for  matrices  of 
Fourier  coefficients,  and  Q  are  the  torus  coordinates  which  are  linear  functions  of 
time,  with  time  derivatives  equal  to  the  basis  frequency  set.  As  noted  before,  Craft 
showed  that  the  optimum  number  of  terms  in  this  Fourier  series  was  between  1500 
and  2000.  Consequently,  1750  terms  were  included  in  the  current  effort. 

The  Fourier  coefficients  can  be  solved  for  by  a  number  of  methods  as  described 
in  [35].  The  method  used  for  this  work  was  the  frequency  cluster  decomposition 
method  studied  and  demonstrated  by  Bordner  [2],  This  method  takes  advantage  of 
the  fact  that  maxima  in  the  power  spectra  occur  at  intervals  equal  to  the  second 
basis  frequency,  cj2,  as  shown  previously  in  Figure  9  and  smaller  peaks  around  these 
are  separated  by  uq.  This  can  be  seen  in  Figure  10.  Consequently,  clusters  of  peaks 
can  be  analyzed  to  determine  their  Fourier  coefficients  simultaneously.  This  method 
simplifies  the  calculations  in  that  it  transforms  the  calculation  of  Fourier  coefficients 
from  solving  a  single  large  linear  system  to  simultaneously  solving  multiple  small 
linear  systems  which  reduces  the  total  amount  of  calculation  required. 

3.6  Transition  to  a  Nearby  Torus 

The  torus  obtained  through  a  first  iteration  of  the  method  described  in  Sec¬ 
tion  3.5  was  a  valid  KAM  torus,  however  it  is  not  the  torus  which  will  describe  the 
motion  of  the  satellite  in  question  because  its  basis  frequencies  were  not  the  desired 
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Partial  Thor  Rocket  Body  Frequency  Power  Spectrum 


Figure  10:  Thor  Rocket  Body  Frequency  Power  Spectrum  Peak  Separation  by  CU3 

frequencies,  those  obtained  via  the  TLE  curve-fits  in  Section  3.2.  This  is  most  likely 
due  to  the  inherent  inaccuracies  of  SGP4,  specifically  the  velocity.  Therefore,  the 
initial  velocity  was  changed  in  such  a  way  that  the  method  of  Section  3.5  yielded  a 
torus  with  the  desired  frequencies.  This  was  done  using  the  following  process. 

The  reference  (first-iteration)  torus  basis  frequencies  are  grouped  into  a  vector, 
lu,  while  the  desired  torus  frequencies  from  the  TLE  data  are  grouped  into  co0.  Then, 
the  frequency  error  can  be  written  as 

did  -> 

A  u  =  uj  —  cu0  =  — ^A  P  (58) 

dP 

where  P  is  a  vector  of  the  reference  torus  momenta  and  A P  is  the  unknown  mo¬ 
mentum  offset  between  the  reference  and  desired  tori.  Equation  58  can  be  solved  for 
A P  which  yields 

AP  =  -A  Acn  (59) 

dP 
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The  Jacobian  matrix  |)=  can  be  found  analytically  using  2BP  orbital  elements  and 
taking  the  torus  momenta  to  be  the  DeLaunay  momenta, 


Pi  =  yfjm 

P-2  =  y rjld\/l  —  e2  cosi  (60) 

P3  =  y/JIdVl  -  e 2 

The  Hamiltonian  can  then  be  written,  through  J2,  as 

K.  _  f  (Pi  -  3 Pi) 

K  2  p2  +  4Pfp| 

which  is  the  same  expression  as  Equation  43.  The  approximate  torus  frequencies, 
again  through  only  J2,  can  be  expressed  as  partial  derivatives  of  the  Hamiltonian 

_  dlC 
oj  =  — ^ 

OP 

and  therefore  the  Jacobian  from  Equation  59  can  be  expressed  as 

du  _  d2IC 
dP  ~  dP2 

and  the  torus  momentum  offsets  can  be  calculated  from  a  known  frequency  offset. 

In  order  to  form  the  new  torus  with  the  desired  basis  frequencies,  a  new  nu¬ 
merical  integration  must  be  carried  out  starting  from  some  new  initial  conditions, 
Xq  =  [qfEo1]7^  Therefore,  the  torus  momentum  offset  must  be  expressed  as  a 
change  in  physical  position  and/or  velocity.  For  this  effort,  the  initial  position  from 
the  TLEs  was  assumed  accurate,  and  therefore  held  constant,  and  the  initial  velocity 
was  allowed  to  change. 
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Let  the  physical  and  torus  state  vectors  be 


*  =  [ 
?=[ 


ft rf 

<r  pt  1 


=  [xy  zvx  vy  vz]T 
=  [MnuP1P2P3}T 


(64) 


Then  the  change  in  the  physical  state  vector  can  be  written  in  terms  of  a  change  in 
the  torus  state  vector  as 


AX 


A  r 
Av 


(65) 


where  AY 


A  QT  A  PT 


T 


.  The  Jacobian 


ax 

dY 


can  be  found  analytically  as 


dX  (dY  dz\  1 
dY  ~  \dZdX ) 


where 

Z  =  [M  VLuj  a  e  i\T 


(66) 


(67) 


is  a  vector  containing  the  classical  orbital  elements.  For  a  detailed  derivation  of 
the  content  of  see  [32],  Then,  setting  the  change  in  initial  position,  A r  =  0, 
Equation  65  can  be  re-written  as  two  linear  equations 


Af  =  0  =  AAQ  +  BAP 


(68) 


Ail  =  CAQ  +  DAP 


where  A,  B,  C,  and  D  are  quadrants  of 


af 

dY 


ax. 

dY' 

A 

C 


B 

D 


(69) 


(70) 
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Given  A P  from  Equation  59,  the  change  in  initial  velocity  needed  to  give  a  desired 
change  in  torus  basis  frequencies  can  be  found  by  first  solving  Equation  68  for  A Q 
and  then  solving  Equation  69  for  An: 

Av  =  (~CA~lB  +  D)  A P  (71) 

Having  now  obtained  a  new  initial  physical  state,  Xq,  the  torus  construction 
process  of  Section  3.5  was  repeated  to  find  a  new  torus  with  updated  basis  frequen¬ 
cies.  These  frequencies  were  again  compared  to  the  desired  TLE  frequencies,  and 
the  frequency  matching  algorithm  was  repeated  until  the  maximum  basis  frequency 
error  was  down  to  O10~12rad/TU . 

3.7  Physical  Coordinate  Extraction 

The  torus  with  frequencies  matching  the  TLE  frequencies  was  then  used  to 
extract  physical  position  and  velocity  as  a  function  of  time.  This  was  done  by 
modifying  a  Matlab  script  initially  developed  by  Capt  Max  Yates  to  read  the  torus 
Fourier  series  hie.  The  original  torus  model  has  coordinates,  Qi,  which  increment 
linearly  in  time  at  rates  equal  to  the  torus  basis  frequencies 


Qi  —  w{t  +  Qio  (72) 

Where  the  Qi0  are  the  value  of  each  coordinate  at  epoch.  However,  as  noted  in 
Section  3.2,  the  TLE  data  showed  that  M,  Q,  and  u  change  as  quadratic  functions  of 
time  and,  since  the  initial  velocity  was  allowed  to  change  in  the  torus  fitting  process, 
the  Qio  differed  slightly  from  the  values  calculated  in  the  TLEs.  The  second-order 
effects  were  taken  into  account  in  the  conversion  from  torus  to  physical  coordinates 
by  calculating  the  torus  coordinates  at  each  time-step  as 


Qi 


0>2  d2  +  W it  +  QiO 


(73) 
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Where  «2i  are  the  2nd-order  coefficients  from  the  TLE  curve  fits.  The  change  in 
Qio  was  addressed  by  manually  changing  the  values  in  the  torus  hie  to  match  the 
a0i  values  in  the  curve  fits.  Having  obtained  the  torus  coordinates,  Equation  57  was 
used  to  transform  to  the  physical  coordinates,  q  =  [. X  Y  Z]T . 

3. 8  Summary 

The  procedure  for  the  current  work  can  then  be  concisely  summarized  in  the 
following  manner: 

1.  Fit  quadratic  curves  to  M,  fi,  and  to  data  from  a  series  of  TLEs.  Set  the  desired 
torus  basis  frequencies  equal  to  the  first-order  curve  fit  coefficients. 

2.  Numerically  integrate  the  orbit  starting  from  an  initial  position  and  velocity 
derived  from  the  Erst  TLE  and  SGP4  algorithms. 

3.  Run  the  torus  construction  algorithm  to  obtain  an  initial  torus  with  basis 
frequencies  near  the  desired  values  from  step  1. 

4.  Perform  an  iteration  of  the  frequency  matching  algorithm  to  match  the  torus 
basis  frequencies  to  the  desired  values. 

5.  Extract  physical  torus  position  data  using  the  torus  Fourier  series.  Calculate 
the  torus  coordinates  within  the  Fourier  series  using  the  quadratic  curve-fits 
from  step  1. 

6.  Compare  the  position  predicted  by  the  torus  to  the  position  at  epoch  of  each 
TLE. 

The  results  of  this  procedure  for  each  test  object  are  presented  and  discussed  in  the 
following  chapter. 
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IV.  Results 


The  process  described  in  Chapter  III  was  attempted  for  each  of  the  four  test  satellites 
with  varying  degrees  of  success.  This  chapter  will  discuss  the  results  from  each  test 
case  and  propose  possible  causes  for  the  difficulties  encountered. 

4-1  Delta  Rocket  Body  No.  1  Results 

4-1.1  TLE  Analysis.  The  TLE  data  for  M,  D,  and  to  can  be  seen  in 

Figures  11  through  13.  The  2nd  order  curves  fit  to  the  data  can  be  seen  in  Table  5. 

The  residuals  for  the  11  and  to  data  are  quite  small,  on  the  order  of  10-4  and  10~3 
respectively,  however  the  residuals  for  the  mean  anomaly  are  between  ±20radians. 
This  poor  fit  is  believed  to  be  due  to  inaccurate  revolution  number  data  in  the  TLE 
set.  This  inaccuracy  is  a  result  of  the  way  the  TLE  data  is  set  up.  Revolution 
number  is  incremented  at  ascending  node  crossing,  and  the  TLE  epoch  time  is  also 
set  to  the  time  when  the  satellite  is  at  the  ascending  node.  Having  epoch  time  and 
revolution  number  tied  to  the  same  point  can  cause  the  revolution  number  to  not  be 
incremented  properly  if  the  position  of  the  satellite  is  off  slightly  from  the  ascending 
node.  Further  analysis  could  be  performed  to  compensate  for  these  errors,  however 
this  analysis  was  not  done  as  part  of  the  current  work. 

Because  of  the  poor  mean  anomaly  curve-fit,  an  attempt  was  not  made  to  fit 
a  torus  to  the  data  for  Delta  Rocket  Body  No.l. 


Table  5:  Delta  Rocket  Body  No.  1  Curve-Fits 


Data 

ao  [rad] 

rad 
ai  TU. 

rad 
cl2  TU2 

M 

6.6068180759878032 

5.36659695481893000 

-1. 30245451 198124E-09 

n 

4.0186766299246193 

-3.03925873725682E-4 

-4.66094156237480E-13 

to 

2.6190404421881577 

3.24842080054474E-4 

-6.36431896081014E-13 

Curve  fit  of  the  form:  X  =  ao  +  a\t  +  02? 
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Mean  Anomaly 
Delta  Rocket  Body  No.1 


Figure  11:  Mean  Anomaly  Data  for  Delta  Rocket  Body  No.l 


Right  Ascension  of  the  Ascending  Node 
Delta  Rocket  Body  No.1 


Residuals  (Data  -  Curve  Fit) 


Figure  12:  RAAN  Data  for  Delta  Rocket  Body  No.l 
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Argument  of  Perigee 
Delta  Rocket  Body  No.1 


Time  (TU) 


Figure  13:  Argument  of  Perigee  Data  for  Delta  Rocket  Body  No.l 


4-2  Delta  Rocket  Body  No.  2  Results 

4-2.1  TLE  Analysis.  The  orbit  of  Delta  Rocket  Body  No. 2  had  an  eccen¬ 
tricity  of  approximately  0.03  and  had  an  apogee  of  only  700km.  The  TLE  data  can 
be  seen  in  Figures  14  through  16.  The  2nd  order  curves  fit  to  the  data  can  be  seen 
in  Table  6. 

The  residuals  for  the  D  and  ui  data  are  again  small,  although  larger  than  those 
for  Delta  Rocket  Body  No.l,  but  the  mean  anomaly  residuals  are  again  the  largest, 
this  time  approximately  ±2  radians.  This  poor  fit  could  be  due  to  the  relatively 
low  altitude  of  the  orbit  (perigee  at  357km)  which  would  cause  this  rocket  body 


Table  6:  Delta  Rocket  Body  No.  2  Curve-Fits 


Data 

ao  [rad] 

rad 
al  TU 

rad 

TTJ2 

M 

7.46732712440206110 

0.88288002293142800 

2.28393350843531E-08 

n 

3.68122129498721940 

-1.10288004166225E-3 

-6.52190507528179E-11 

U ) 

0.17914312485328654 

1.88174783746889E-3 

1.11759586936959E-10 

Curve  fit  of  the  form:  X  =  ao  +  ait  +  a2i2 
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Mean  Anomaly 
Delta  Rocket  Body  No.2 


Time  (TU) 


Figure  14:  Mean  Anomaly  Data  for  Delta  Rocket  Body  No.2 


Right  Ascension  of  the  Ascending  Node 
Delta  Rocket  Body  No.2 


Time  (TU) 


Figure  15:  RAAN  Data  for  Delta  Rocket  Body  No.2 
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120 


Argument  of  Perigee 
Delta  Rocket  Body  No.2 


Time  (TU)  x104 

Figure  16:  Argument  of  Perigee  Data  for  Delta  Rocket  Body  No.2 

to  encounter  significant  air  drag.  And,  while  the  difference  between  apogee  and 
perigee  is  only  approximately  350km,  because  the  effect  of  air  drag  increases  roughly 
exponentially  with  decreasing  altitude  there  would  still  be  a  significant  difference  in 
the  amount  of  air  drag  felt  over  the  course  of  an  orbit. 

Again,  because  of  the  poor  curve  fits,  a  torus  was  not  fit  to  the  data  for  Delta 
Rocket  Body  No.2. 

4-3  Hubble  Space  Telescope  Results 

4-3.1  TLE  Analysis.  The  Hubble  Space  Telescope  is  in  a  nearly  circular  or¬ 
bit  with  eccentricity  of  approximately  0.0004  at  an  altitude  of  approximately  560knr. 
In  addition,  while  the  Rocket  Bodies  analyzed  are  large  hollow  cylinders  with  a  large 
cross-sectional  area  to  mass  ratios  (^),  the  HST  is  filled  with  optics,  cameras,  bat¬ 
teries  and  other  equipment  which  give  it  much  smaller  j~r  Because  of  these  factors, 
it  was  expected  the  air-drag  effects  would  be  less  pronounced  and  nearly  uniform 
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throughout  the  orbit.  TLE  data  can  be  seen  in  Figures  17  through  19  while  the  2nd 
order  curves  fit  to  the  data  can  be  seen  in  Table  7. 

Table  7:  Hubble  Space  Telescope  Curve-Fits 


Data 

a0  [rad] 

rad 
ai  xu 

'  rad 
d2  TU2 

M 

3.48798565697538270 

0.88019778693820000 

1.01689070044144E-09 

n 

0.71456401402421699 

-1.06265838887821E-3 

-2.81526562126039E-12 

uj 

2.71458280469347900 

1.73080477457464E-3 

-3.66447396216027E-12 

Curve  fit  of  the  form:  X  =  qq  +  a±t  + 


xio' 


Mean  Anomaly 
Hubble  Space  Telescope 


Time  (TU) 


Figure  17:  Mean  Anomaly  Data  for  the  HST 


This  time,  residuals  for  are  still  small,  but  both  mean  anomaly  and  u  have 
residuals  of  approximately  ±0.1  radian.  These  residuals  may  again  be  influenced  by 
the  randomness  of  air  drag  caused  by  changes  in  the  earth’s  atmosphere.  While  the 
residuals  may  seem  small,  it  is  important  to  remember  that  an  error  of  0.1  radians 
in  mean  anomaly  translates  to  a  position  difference  of  nearly  700km  at  the  HST’s 
orbital  altitude.  Because  of  this,  a  torus  accurate  to  even  10s  of  kilometers  in  position 
would  not  be  possible  with  the  current  method,  however  the  torus  fitting  process  was 
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Right  Ascension  of  the  Ascending  Node 
Hubble  Space  Telescope 


Residuals  (Data  -  Curve  Fit) 


Figure  18:  RAAN  Data  for  the  HST 


still  attempted  to  determine  if  it  would  be  accurate  to  expected  precision  based  on 
the  residuals. 
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Argument  of  Perigee 
Hubble  Space  Telescope 


Time  (TU)  x104 

Figure  19:  Argument  of  Perigee  Data  for  the  HST 

4-3.2  Torus  Fitting.  The  desired  frequencies  for  the  HST  torus  were  taken 
from  the  first  order  terms  in  the  curve  fits  and  can  be  seen  in  Table  8. 


Table  8:  Hubble  Space  Telescope  Desired  Torus  Frequencies 


UJl 

rad 

.TU. 

OJ  2 

rad 

.TU. 

u  3 

rad 

.TU. 

8.80I977869382E-1 

-5.98963I4969997IE-2 

I.73080477457464E-3 

The  torus  construction  process  of  Section  3.5  was  attempted,  but  failed  on 
the  first  attempt.  The  algorithm  to  identify  the  basis  frequencies  in  the  numerically 
integrated  data  failed  to  correctly  identify  cu3  starting  from  the  J2  estimate.  In  an 
attempt  to  solve  this  problem,  the  frequency  power  spectrum  was  analyzed  manually 
to  get  a  better  estimate  of  cu3.  Starting  from  this  estimate,  the  frequency  basis  set 
improved,  but  still  had  residuals  on  the  order  of  lOA— 5  which  showed  that  the  correct 
basis  frequencies  had  still  not  been  identified.  It  is  believed  that  the  cause  of  this 
difficulty  may  be  the  small  eccentricity  of  the  HST’s  orbit.  I11  classical  perturbation 
theory,  small  eccentricities  lead  to  singularities  in  the  series  expansions  of  the  orbital 
elements.  These  singularities  are  dealt  with  in  the  SGP4  algorithms  by  discarding 
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the  terms  that  cause  problems  [12].  The  torus  construction  method,  however,  does 
not  neglect  these  terms  and  the  low  eccentricity  of  the  HST’s  orbit  seems  to  be  below 
the  useful  tolerance  of  the  current  method. 

4-4  Thor  Rocket  Body  Results 

4-4-1  TLE  Analysis.  The  Thor  Rocket  Body  is  in  an  orbit  with  an  eccen¬ 
tricity  of  approximately  0.02  with  apogee  at  an  altitude  of  1203  km  and  perigee  at 
921  km.  With  this  orbit,  it  was  expected  that  air  drag  effects  would  be  less  pro¬ 
nounced  than  what  was  seen  with  Delta  rocket  body  No. 2  or  the  HST.  The  TLE 
data  can  be  seen  in  Figures  20  through  22  and  the  curve  fits  are  shown  in  Table  9. 


Mean  Anomaly 
Thor  Rocket  Body 


Time  (TU) 


Figure  20:  Mean  Anomaly  Data  for  the  Thor  Rocket  Body 

This  time,  all  residuals  were  small  with  a  maximum  value  of  approximately  5 E  —  3 
seen  in  the  mean  anomaly  residuals.  Because  of  this,  the  torus  fitting  process  was 
attempted. 
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Right  Ascension  of  the  Ascending  Node 
Thor  Rocket  Body 


Time  (TU)  x104 

Figure  21:  RAAN  Data  for  the  Thor  Rocket  Body 


Table  9:  Thor  Rocket  Body  Curve-Fits 


Data 

a0  [rad] 

rad 
ai  TU 

'  rad 
cA2  TU2 

M 

4.8169895251603805 

0.79372722076537100 

5.93349302063384E-11 

n 

1.9555018406270159 

-8.36935587639413E-4 

-1.45394333352382E-13 

IV 

1.5043570873484007 

1.36750321345161E-3 

4.85354594310724E-13 

Curve  fit  of  the  form:  X  =  a®  +  a\t  +  a-fi1 
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Argument  of  Perigee 
Thor  Rocket  Body 


Time  (TU) 


xIO4 


Figure  22:  Argument  of  Perigee  Data  for  the  Thor  Rocket  Body 
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4-4-2  Torus  Fitting.  Based  on  the  curve  fits,  the  desired  torus  basis  fre¬ 
quencies  were  determined  and  are  shown  in  Table  10.  Starting  from  the  initial  po- 


Table  10:  Thor  Rocket  Body  Desired  Torus  Frequencies 


Ui 

'rad' 

lTU\ 

U)2 

'rad' 

LTUl 

UJs 

'rad' 

YtUI 

7.9372722076537E-1 

-5.96705921687582E-2 

1.3675032134516E-3 

sition  and  velocity  determined  using  the  SGP4  code  and  the  first  TLE,  the  process 
of  creating  a  torus  and  matching  the  basis  frequencies  to  the  desired  set  described 
in  Sections  3.5  and  3.6  was  accomplished.  A  total  of  four  iterations  of  the  frequency 
matching  algorithm  were  necessary  to  match  the  torus  basis  frequencies  to  within  the 
desired  tolerance.  Table  11  shows  the  progression  of  torus  frequencies  and  Table  12 
shows  the  initial  velocity  changes  made  at  each  iteration  of  the  frequency  matching 
process. 


Table  11:  Torus  Frequency  Matching  for  Thor  Rocket  Body 


, ,  'rad' 
U 1  TU 

,  ,  'rad' 

UJ2  TU 

.  ,  'rad' 
UJ3  TU 

Initial  Torus 

7.93732038441640E-01 

-5.96704909422124E-02 

1.36749993586527E-03 

Iteration  1  Torus 

7.93718780146312E-01 

-5.96705786750190E-02 

1.36747666904746E-03 

Iteration  2  Torus 

7.93727192361304E-01 

-5.96705924650170E-02 

1.36750386446849E-03 

Iteration  3  Torus 

7.93727221037809E-01 

-5.96705921734135E-02 

1.36750322508483E-03 

Iteration  4  Torus 

7.93727220771728E-01 

-5.96705921687751E-02 

1.36750321348700E-03 

u i  residual 

'rad' 

Ytu\ 

0J2  residual 

'rad' 

Ytu\ 

0J3  residual 

rad' 

Ytu\ 

Iteration  1  Torus 

8.44061905791449E-06 

-1.34937392023970E-08 

2.65444041400165E-08 

Iteration  2  Torus 

2.84040659881413E-08 

2.96258795273729E-10 

-6.51016889919101E-10 

Iteration  3  Torus 

-2.72439071302699E-10 

4.65529698123746E-12 

-1.16332299696098E-11 

Iteration  4  Torus 

-6.35802521742335E-12 

1.68962066560141E-14 

-3.53998514529552E-14 

Table  12:  Changes  to  Initial  Velocity  in  Thor  Rocket  Body  Torus  Fitting  Process 


V 

v  x  TU 

V 

v  y  tu 

V 

v  z  Tu 

Initial  Velocity 

-9.98015949328921E-2 

-8.11259141324669E-1 

4.39329318458694E-1 

Iteration  1  AV 

2.45940386540720E-3 

-2.58556420032001E-4 

8.51464239486448E-5 

Iteration  2  AV 

-1.09486652807646E-4 

1.50534148714208E-5 

-3.35944758461485E-6 

Iteration  3  AV 

-3.85930049076970E-6 

5.54734206659374E-7 

1.45063141307529E-7 

Iteration  4  AV 

-3.20433729724457E-8 

5.37472590607683E-9 

3. 03779212051 191E-9 

Total  AV 

-2.34602586873582E-3 

2.42942896228060E-4 

-8.19350772974792E-5 
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The  magnitude  of  the  total  velocity  change  was  approximately  0.25  percent  of 
the  initial  velocity  magnitude  obtained  from  the  TLEs  and  SGP4.  This  corresponds 
to  approximately  19  m/s  and  is  well  within  the  uncertainty  in  velocity  predictions 
using  SGP4. 


4-4-3  Position  Comparison.  After  a  torus  was  formed  with  the  desired 
basis  frequency  set,  the  torus  coordinates  were  transformed  to  physical  coordinates 
using  Equation  57  and  a  Matlab  script.  This  physical  position  was  then  compared 
to  the  numerically  integrated  position  to  determine  if  the  torus  matched  the  data  it 
was  fit  to.  The  results  of  this  comparison  can  be  seen  in  Figure  23. 


Figure  23:  Comparison  of  Numerically  Integrated  Position  to  Torus  Position,  1st 
order  Qi 


The  torus  position  matches  the  numerically  integrated  position  to  within  ap¬ 
proximately  60  meters  in  any  of  the  three  axes  over  6  months.  The  abrupt  change 
in  the  error  growth  rate  seen  at  approximately  170  days  is  the  result  of  the  torus 
coordinate  calculations  reaching  the  limit  of  double  precision  accuracy. 


57 


The  position  predicted  by  the  torus  was  then  calculated  at  each  TLE  epoch 
time  and  compared  to  the  SGP4  position.  This  was  initially  done  using  only  the 
first  order  curve-fit  terms  (the  basis  frequencies)  to  update  the  torus  coordinates, 
Qi,  and  the  initial  torus  coordinate  values,  Qm,  calculated  in  the  torus  construction 
process.  The  results  of  this  comparison  can  be  seen  in  Figures  24  and  25. 


Residuals:  SGP4  Position  -  Torus  Position  (1st  Order) 


Figure  24:  Residuals: 


SGP4  Position  -  Torus  Position,  1st  order  Qi 
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Magnitude  of  Position  Difference:  SGP4  Position  -  Torus  Position  (1st  Order) 


Figure  25:  Error  Magnitude:  SGP4  Position  -  Torus  Position,  1st  order  Qi 

It  is  evident  that  the  torus  in  this  configuration  provides  a  poor  position  pre¬ 
diction  with  the  error  in  position  going  to  nearly  1500km  after  18  months.  This, 
however,  is  expected  as  the  TLE  analysis  showed  that  the  torus  coordinates,  M,  hi, 
and  uj  increase  as  quadratic,  2nd  order  functions  of  time,  not  first  order.  To  com¬ 
pensate  for  this,  the  2nd  order  curve-fit  coefficients  were  included  in  the  calculation 
of  the  Q's  at  each  time-step  and  the  torus  position  was  again  compared  to  the  SGP4 
position.  The  results  of  this  comparison  can  be  seen  in  Figures  26  and  27. 
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Residuals:  SGP4  Position  -  Torus  Position  (2nd  Order) 
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6:  Residuals:  SGP4  Position  -  Torus  Position,  2nd  order  Qi 


Magnitude  of  Position  Difference:  SGP4  Position  -  Torus  Position  (2nd  Order) 
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Figure  27:  Error  Magnitude:  SGP4  Position  -  Torus  Position,  2nd  order  Q 


These  plots  show  that  adding  the  2nd  order  terms  to  the  torus  coordinates 
drastically  improve  the  torus  position  predictions,  with  the  maximum  error  now 
being  approximately  110km,  or  less  than  10  percent  of  the  previous  error.  It  is 
also  important  to  note  that  the  linear  portion  of  the  error  growth  has  been  nearly 
eliminated  leaving  only  periodic  oscillations. 

In  an  effort  to  further  refine  the  torus  position  prediction,  the  a o  terms  for  the 
curve-fits  were  used  to  replace  the  Qio  values  calculated  in  the  torus  construction 
process.  This  replacement  is  shown  in  Equation  74. 

Qio  =  Mo 

Q20  =  O0  -  GMST0  (74) 

Q  30  — 

The  torus  position  was  again  compared  to  the  SGP4  position  and  can  be  seen  in 
Figures  28  and  29. 


Residuals:  SGP4  Position  -  Torus  Position  (2nd  Order) 


Figure  28:  Residuals:  SGP4  Position  -  Torus  Position,  2nd  order  Qj,  Q0i  from  curve 
fits 
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Magnitude  of  Position  Difference:  SGP4  Position  -  Torus  Position  (2nd  Order) 
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Figure  29:  Error  Magnitude:  SGP4  Position  -  Torus  Position,  2nd  order  Qi , 
from  curve  fits 

This  comparison  shows  a  number  of  notable  results.  First,  the  maximum  error 
is  now  down  to  approximately  78km;  a  29  percent  reduction  from  the  error  obtained 
using  the  initial  Qi o  values.  Second,  there  is  now  error  present  at  t=0.  This  initial 
error  is  due  to  the  fact  that  the  values  of  the  curve  fits  for  M,  12,  and  to  at  t=0  are 
not  exactly  equal  to  the  values  in  the  first  TLE.  Finally,  the  amplitude  of  oscillation 
in  the  X,  Y,  and  Z  error  has  decreased.  This  is  shown  in  Table  13. 


Table  13:  Comparison  of  Position  Error  Oscillation  Amplitude 


X  Error  Oscillation 
Amplitude  [km] 

Y  Error  Oscillation 
Amplitude  [km] 

Z  Error  Oscillation 
Amplitude  [km] 

Initial  Q^ 

155  -  180 

85  -  120 

35  -  42 

Curve  fit  Qi o 

85  -  125 

55  -  86 

37-39 

While  this  representation  of  the  torus  does  drastically  reduce  the  error  in  posi¬ 
tion  from  the  initial  torus  with  first-order  Q’s,  it  still  does  not  provide  anywhere  near 
the  accuracy  that  would  be  required  for  operational  use.  In  an  effort  to  determine 
the  cause  of  the  error,  the  periodic  nature  of  the  error  was  examined. 
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Figure  28  shows  that  the  X  and  Y  error  oscillations  are  nearly  in-phase,  while 
the  Z  oscillations  are  offset  by  180  degrees.  It  is  also  interesting  to  note  that  the  Z 
error  oscillations  seem  to  have  only  one  frequency,  while  the  X  and  Y  oscillations 
have  at  least  two  distinct  frequencies.  The  period  between  peaks  in  the  Z  data  is 
approximately  42.2  to  43.3  days,  which  is  also  the  period  between  the  major  peaks 
in  the  X  and  Y  data.  This  period  corresponds  to  the  period  of  0J3,  which  is  42.905 
days.  To  further  verify  this,  a  FFT  was  performed  on  the  position  error  data  to 
analyze  its  frequency  content.  Figure  30  shows  the  frequency  power  spectrum  over 
the  entire  range  of  frequencies  studied. 


Figure  30:  Frequency  Power  Spectrum  for  Torus  Position  Error 

The  spectrum  is  quite  noisy,  but  clearly  has  maximum  power  in  the  smaller 
frequencies.  The  two  points  in  the  spectrum  where  the  oscillations  appear  to  tem¬ 
porarily  damp  out  are  centered  around  frequencies  equal  to  the  torus’  uq  and  2uq. 
Figure  31  shows  detail  around  uj\.  It  is  clear  that  there  is  no  peak  at  cui,  but  there 
are  small  peak  clusters  at  intervals  of  ±a;2  separated  by  oj3  in  the  X  and  Y  spectra. 
Figure  32  shows  detail  around  the  torus’  co3  value  and  shows  that,  as  expected,  all 
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three  coordinates  have  maximum  power  at  this  frequency.  Also,  as  expected,  the  Z 
coordinate  power  spectrum  drops  off  sharply  around  the  CU3  value  while  the  X  and 
Y  spectra  have  other  maxima  of  similar  magnitude  nearby. 


Figure  31:  Frequency  Power  Spectrum  for  Torus  Position  Error  -  ui\  detail 


64 


Frequency  Spectrum  of  Torus  Position  Error  -  Detail  around  w. 


Figure  32:  Frequency  Power  Spectrum  for  Torus  Position  Error  -  uq  detail 

Figure  33  shows  detail  around  the  torus’  002  value  and  clearly  shows  three  more 
peaks  in  X  and  Y  at  positions  corresponding  to  CO2  and  0J2  uj:i .  The  Z  coordinate  is 
not  expected  to  have  peaks  at  these  locations  because  0J2,  or  Q,  has  minimal  impact 
on  Z-position.  (Recall  that  the  numerically  integrated  Z-Coordinate  data  did  not 
have  peaks  at  [1  1  1]  or  [1  -1  1]). 
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Frequency  Spectrum  of  Torus  Position  Error  -  Detail  around  w2 


Figure  33:  Frequency  Power  Spectrum  for  Torus  Position  Error  -  u2  detail 

This  analysis  suggests  that  the  torus  construction  algorithm  does  not  ade¬ 
quately  characterize  the  impacts  of  cu3  and,  to  a  smaller  extent,  u2  on  the  satellite’s 
position.  In  order  to  determine  if  this  deficiency  was  caused  by  too  few  terms  in  the 
Fourier  series,  another  torus  was  constructed  containing  additional  terms,  how¬ 
ever  the  position  error  measured  against  the  TLE  and  SGP4  data  was  unchanged. 
Because  of  this,  it  is  believed  that  the  algorithm  used  to  calculate  the  Fourier  series 
coefficients  through  simultaneously  analyzing  clusters  of  peaks  is  not  working  prop¬ 
erly.  In  this  method,  clusters  of  peaks  are  studied  in  order  of  decreasing  magnitude. 
Each  peak  cluster  is  analyzed  in  order  to  determine  its  Fourier  coefficients  and  sub¬ 
sequently  removed  from  the  data.  This  technique  is  used  in  the  hope  that  removing 
larger  peak  clusters  will  allow  more  precise  analysis  of  smaller  clusters  which  would 
otherwise  be  obscured  by  their  larger  neighbors.  The  holes  in  the  power  spectrum  of 
the  position  error  at  oj\  and  2uj\  demonstrate  that  the  larger  oji  peak  clusters  have 
been  successfully  eliminated,  however  the  appearance  of  peaks  surrounding  u>2  and 
cu3  suggest  that  these  smaller  peak  clusters  have  not  been  removed  properly  from 
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the  data.  In  his  research,  Bordner  also  experienced  problems  in  properly  charac¬ 
terizing  the  impact  of  u>3  [2],  He  concluded  that  these  problems  stemmed  from  two 
sources.  First,  that  if  the  orbital  data  analyzed  spanned  too  short  a  time  period 
then  the  effect  of  u3  could  not  be  captured  due  to  it’s  long  period,  and  second  that 
low  eccentricity  orbits  would  also  encounter  difficulty  clue  to  small  divisors.  It  was 
believed  that  the  current  work  would  avoid  each  of  these  issues  as  the  time  period  of 
data  used  to  construct  the  torus  is  approximately  13  times  the  period  of  cv3  and  the 
eccentricity  of  the  orbit  is  approximately  0.02.  Further  research  is  needed  to  verify 
that  the  current  difficulty  is  another  manifestation  of  problems  faced  by  Bordner,  or 
if  it  is  a  new  issue  entirely. 
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V.  Conclusions 


Chapter  IV  shows  a  number  of  promising  results.  First,  the  HST  and  Thor  Rocket 
Body  test  cases  show  that  it  is  possible  to  extract  accurate  KAM  torus  basis  frequen¬ 
cies  from  observational  data  in  the  form  of  TLEs.  Second,  a  method  of  making  slight 
changes  to  torus  basis  frequencies  was  demonstrated  by  equating  a  desired  change 
in  basis  frequencies  to  a  small  change  in  initial  velocity.  This  capability  allows  the 
current  torus  construction  method  to  be  used  in  tandem  with  the  TLE  frequency 
identification  routine  to  obtain  a  torus  with  the  desired  basis  frequencies  for  a  given 
satellite. 

The  difficulties  encountered  with  TLE  curve  fitting  for  the  two  Delta  Rocket 
Bodies  show  that  there  may  be  a  only  a  limited  subset  of  orbital  regimes  for  which 
basis  frequency  extraction  from  TLEs  can  be  useful;  Specifically,  orbital  regimes  with 
moderate  inclination  and  eccentricity  where  air  drag  is  minimal  and  nearly  constant 
throughout  the  orbit.  The  difficulty  seen  in  the  torus  construction  process  for  the 
HST  and  the  poor  characterization  of  and  CU3  in  the  Thor  Rocket  Body  test  case 
reinforces  the  fact  that  small  eccentricity  orbits  cannot  be  modeled  to  the  desired 
degree  of  accuracy  with  the  current  KAM  torus  construction  method. 

5.1  Recommendations  for  Further  Study 

First,  an  investigation  into  the  failure  of  the  peak-cluster  decomposition  method 
of  calculating  Fourier  Coefficients  seen  in  the  Thor  Rocket  Body  test  case  is  needed. 
This  test  case  should  be  continued,  either  using  an  updated  peak-cluster  method  or 
using  the  original  Laskar  method  of  analyzing  each  peak  individually.  This  investi¬ 
gation  should  be  completed  to  show  the  true  accuracy  of  the  torus  model. 

Next,  a  broader  survey  is  needed  of  satellites  in  a  variety  of  orbits.  This  analysis 
should  be  conducted  to  show  the  limits  of  eccentricity,  inclination,  and  orbital  period 
that  form  the  boundaries  of  the  ‘orbital  box’  in  which  the  extraction  of  torus  basis 
frequencies  from  TLEs  is  possible.  Another  study  could  be  done  to  determine  the 


time  period  of  TLE  data  that  is  necessary  to  obtain  accurate  basis  frequencies,  ft 
would  also  be  beneficial  to  investigate  the  possibility  of  reformulating  the  current 
torus  construction  method  in  Poincare  elements.  If  successful,  this  could  eliminate 
the  problems  currently  faced  in  applying  the  KAM  torus  model  to  earth  satellites  in 
low  eccentricity  orbits. 

This  research  on  the  contribution  of  TLE  data  to  KAM  torus  construction 
should  also  be  incorporated  with  the  research  currently  being  done  by  Capt’s  Hagen 
and  Yates.  Their  projects  will  add  air  drag  and  the  effect  of  the  moon  to  the  torus 
model  as  well  as  a  method  of  moving  between  nearby  tori  with  slightly  different 
momenta.  This  integration  could  increase  the  accuracy  of  torus  position  prediction 
and  possibly  allow  the  construction  of  accurate  tori  for  satellites  with  larger  residuals 
in  their  TLE  curve-fits  due  to  their  orbits  or  periodic  station-keeping  maneuvers.  A 
study  could  also  be  conducted  to  determine  how  much  of  the  residuals  seen  in  the 
TLE  curve-fits  are  due  to  stochastics,  such  as  variation  in  air-drag,  and  how  much 
is  due  to  other  perturbations  that  have  not,  or  cannot,  be  modeled  easily. 

Finally,  an  analysis  should  be  conducted  to  compare  the  accuracy  of  a  KAM 
torus  to  the  SGP4  propagation  algorithms.  This  analysis  could  show  the  possible 
benefits  of  using  a  KAM  torus  model  in  place  of  the  current  SGP4  model  in  predicting 
future  satellite  positions. 
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